Subroutine SIMULATION(EXP_L_prime,B, RET_PRIME, jExperimentNumber )

Use PARAMETERS
USE Functions_libs
USE Interpolate_func
USE random

Implicit NONE

Integer :: jzt_int, jo_lagged, juntreated, joption_pay, jr_rep, jtype2_new, jdp_new, jdu_new, js_new,  I_pay,  jhc_ret, jmc, jtrpay, I_treat, jos, jtype4, jtype4_p, jno, jBenchmark_exp, jh_alt, jr_alt, jr_p_alt, jexp_new, jMEcat_p, jMEcat, jExperimentNumber , jspouse,counter, jzhc_p,jzhc, jmar,jmar_p, jtype3, jtype3_p, jzt, jzt_p,jsample,jas,jas_p,jsav,je,jh,jh_p, jr,jr_p,ji, jage,j, jht, jyears, jfix, jtype1,jtype2,jtype2_p,js,jdp,jdu,jx, jo,jo_p, jins, jhrs , jx_age, jo_act , jsav_act, kk  , jsav_upper , TR_indic, jinv_upper  , jsav_upper_ret  , jhk, jinc

Real (Kind=dd), dimension (n_as,n_type1,n_type3,n_hc_grid,n_ret-1,n_o,n_mar,2)::  EXP_L_prime  !value functions    
Real (Kind=dd), dimension (n_o,n_type4,n_zt,n_as,n_type3,n_type1,n_hc_grid,n_ret-1)::  B !value function while working. (1=not have option to treat and not pay, 2=have option)
Real (Kind=dd), dimension (n_h,n_r,n_o,n_mar,n_zhc)::  EXP_L_prime_INT

Real (Kind=dd), dimension (n_hc_ret,n_MEcat,n_as,n_type1,n_h,n_r,n_type2,n_mar):: RET !value functions when retired
Real (Kind=dd), dimension (n_hc_ret,n_MEcat,n_as,n_type1,n_h,n_r,n_type2,n_ret : n_age,n_mar):: RET_PRIME !value functions when retired
Real (Kind=dd), dimension (n_MEcat,n_h,n_r,n_type2,n_mar):: RET_PRIME_INT


Real  (Kind=dd) :: Resources, prob_s, prob_d, income_bt, income, base_1, base_2 , Total_tax, actual_earn, tax_SS, tax_med, V_death, V_c_all_trpay, V_cf_trpay
Real (Kind=dd) :: last_value, maxB, valret   , check
Real (Kind=dd) :: UTI,x, EXPECTED_RET_PRIME 

! these following things depend on experiment. For example, in experiment 12, no one gets a d shock at 60 (36). We want to store the human capital and zt shocks. so we store them for the start of ages 61, 62, 63 and 64 (4 of them). For ages that are younger, we need to store more, up to end of working life. 
Real (Kind=dd), dimension (1:34)::  HC_hypo_9 , zt_hypo_9, o_hypo_9 ! this is for reading the hypothetical HC for each person, as if they did not have a serious shock.
Real (Kind=dd), dimension (1:24)::  HC_hypo_10 , zt_hypo_10, o_hypo_10 ! this is for reading the hypothetical HC for each person, as if they did not have a serious shock.
Real (Kind=dd), dimension (1:14)::  HC_hypo_11 , zt_hypo_11 , o_hypo_11! this is for reading the hypothetical HC for each person, as if they did not have a serious shock.
Real (Kind=dd), dimension (1:4)::   HC_hypo_12 , zt_hypo_12, o_hypo_12 ! this is for reading the hypothetical HC for each person, as if they did not have a serious shock.

Real (Kind=dd), dimension (1:443)::  ran_hypo ! these are the saved random draws from the simulation without shocks. - fill in those with shocks. 6+11*39+5 (6 extra for age 1, and at age 40, we only need to draw 5 shocks that are relevant for decisions at age 40 (64). All other shocks are not relevant until age 65).

Real (Kind=dd):: CON, TR, ASP , MAX_SAV  , as_opt, TR_opt ,   INT_sav_upper_ret   , INT_sav_upper, EXPECTED_B_PRIME
Real (Kind=dd), dimension (n_o) :: EXPECTED_B


Real (Kind=dd) :: jx_int_p, prob_sim
Real (Kind=dd), dimension (n_zhc) :: HK_INTERPOL_sim

Real (Kind=dd):: Z_INT_ALT , Z_INT_OFFER_ALT , temp,  jcons_act, jas_p_act, jtr_act  ,jassets  , jx_int, Z_INT , Z_INT_OFFER , B_act , RET_act  , jx_age_interpol  , EXPECTED_RET_PRIME_act, EXPECTED_B_PRIME_act
Integer :: bb,d,o,f   , gr, gg
Real (Kind=dd) :: rmaxipzt  , rmaxipzhc
Real (Kind=dd), dimension (2,n_MEcat,n_type2,n_o) ::  B_int

Real (Kind=dd), dimension(3)::  CON_sim,  ASP_sim, B_sim,   TR_sim
Real (Kind=dd), dimension(2):: income_pay, Total_tax_pay, base_2_pay, base_1_pay
                    
Integer :: indicator_ret, idum, ind_death, intepsi
Real (Kind=dd), dimension(n_mar) :: rmaxiTPzp
Real (Kind=dd) ::  rmaxiIDzp
Real (Kind=dd) :: rx, prob, ry, rb, rd, rf, ri , ro, rinv, i
Real (Kind=dd) :: epsi,weight_epsi, zp_step


1008 Format(I5,1x,I2, 1x,I1,1x,I1,1x,I1,1x,I1,1x,I1,1x,I1,1x,I1,1x,I1,1x,I1,1x,I1,1x,I1,1x,F6.2 , 1x, F6.2 , 1x, F9.2 , 1x,F9.2, 1x,F9.2, 1x,F10.2, 1x,F10.2, 1x, F9.2, 1x, F8.2 , 1x, F10.2, 1x, F4.1, 1x,F9.2, 1x, F10.2,1x,I1,1x, F9.2,1x,F10.2,1x,F10.2,1x,F10.2,1x,F10.2, 1x, I1, 1x, I1, 1x, I1, 1x, I1, 1x, I1 , 1x, F10.6)
      ! file for writing human capital  - write for all ages from 30 onwards, 40, 50, 60 if no shocks at these ages and no deterioration in health
1004  Format(F16.11)  ! this is for human capital
1005 Format(F16.11) ! this is for random shocks 
1006 Format(I1)  ! zt shocks are 1, 2, or 3
     
     
!Grid for savings is finer in the simulation     
sav_sim(:)=0.0d0
open(258, file='External_Parameters/savings_grid_sim.txt')
do jsav=1, n_sav_sim
   read (258,*) sav_sim(jsav)
end do 
close (258)
sav_sim(:)=sav_sim(:)/5200.0d0     
  

! sav_grid: useful for interpolation to find the maximum savings
sav_grid_sim(:)=0.0d0
 do jsav=1, n_sav_sim
 sav_grid_sim(jsav) = dble(jsav)
 end do
     
! NUMBER OF PEOPLE OF EACH TYPE SIMULATED
do jtype1=1, n_type1
n_samp(jtype1)= dble(n_sample)*In_type1(jtype1)  
end do


!***********************************************************************
! Simulating the economy
!***********************************************************************

! jx_age is an integer number on a constructed grid for human capital and age, e.g., 1 (age=1, hk=1), 2 (age=2, hk=1), 3 (age=2, hk=2) etc
! jx_age_interpol is a real number that can take values in between those of jx_age.  E.g., 4.5 (age=3, hk=1.5, meaning .5 yrs of full time work)
  
 do jBenchmark_exp= 29, 32 ! we need experiments 5, 6, 9-16, 18, 27, 29-32   CHANGE THIS AS NEEDED
     open(14, file='Simulated_Data/Random_numbers0.txt')  !random numbers saved

if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.0 ) then
open(12, file='Simulated_Data/Simulated_data0.txt')
else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.1 ) then
open(12, file='Simulated_Data/Simulated_data01.txt')
else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.2 ) then
open(12, file='Simulated_Data/Simulated_data02.txt')
else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.3 ) then
open(12, file='Simulated_Data/Simulated_data03.txt')
else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.4 ) then
open(12, file='Simulated_Data/Simulated_data04.txt')
else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.5 ) then
open(12, file='Simulated_Data/Simulated_data05.txt')
else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.6 ) then
open(12, file='Simulated_Data/Simulated_data06.txt')
else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.7 ) then
open(12, file='Simulated_Data/Simulated_data07.txt')
else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.8 ) then
open(12, file='Simulated_Data/Simulated_data08.txt')
else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.9 ) then
open(12, file='Simulated_Data/Simulated_data09.txt')
open(13, file='Simulated_Data/Simulated_HC09.txt')
open(15, file='Simulated_Data/Random_jzt09.txt')
open(16, file='Simulated_Data/Simulated_offers09.txt')
else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.10 ) then
open(12, file='Simulated_Data/Simulated_data010.txt')
open(13, file='Simulated_Data/Simulated_HC010.txt')
open(15, file='Simulated_Data/Random_jzt010.txt')
open(16, file='Simulated_Data/Simulated_offers010.txt')
else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.11 ) then
open(12, file='Simulated_Data/Simulated_data011.txt')
open(13, file='Simulated_Data/Simulated_HC011.txt')
open(15, file='Simulated_Data/Random_jzt011.txt')
open(16, file='Simulated_Data/Simulated_offers011.txt')
else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.12 ) then
open(12, file='Simulated_Data/Simulated_data012.txt')
open(13, file='Simulated_Data/Simulated_HC012.txt')
open(15, file='Simulated_Data/Random_jzt012.txt')
open(16, file='Simulated_Data/Simulated_offers012.txt')
else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.13 ) then
open(12, file='Simulated_Data/Simulated_data013.txt')
else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.14 ) then
open(12, file='Simulated_Data/Simulated_data014.txt')
else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.15 ) then
open(12, file='Simulated_Data/Simulated_data015.txt')
else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.16 ) then
open(12, file='Simulated_Data/Simulated_data016.txt')
else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.17 ) then
open(12, file='Simulated_Data/Simulated_data017.txt')
open(13, file='Simulated_Data/Simulated_HC09.txt')
open(15, file='Simulated_Data/Random_jzt09.txt')
else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.18 ) then
open(12, file='Simulated_Data/Simulated_data018.txt')
open(13, file='Simulated_Data/Simulated_HC010.txt')
open(15, file='Simulated_Data/Random_jzt010.txt')
else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.19 ) then
open(12, file='Simulated_Data/Simulated_data019.txt')
open(13, file='Simulated_Data/Simulated_HC011.txt')
open(15, file='Simulated_Data/Random_jzt011.txt')
else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.20 ) then
open(12, file='Simulated_Data/Simulated_data020.txt')
open(13, file='Simulated_Data/Simulated_HC012.txt')
open(15, file='Simulated_Data/Random_jzt012.txt')
else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.27 ) then
open(12, file='Simulated_Data/Simulated_data027.txt')
open(16, file='Simulated_Data/Simulated_offers010.txt')

else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.29 ) then
open(12, file='Simulated_Data/Simulated_data029.txt')
open(13, file='Simulated_Data/Simulated_HC09.txt')
open(15, file='Simulated_Data/Random_jzt09.txt')
open(16, file='Simulated_Data/Simulated_offers09.txt')

else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.30 ) then
open(12, file='Simulated_Data/Simulated_data030.txt')
open(13, file='Simulated_Data/Simulated_HC010.txt')
open(15, file='Simulated_Data/Random_jzt010.txt')
open(16, file='Simulated_Data/Simulated_offers010.txt')

else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.31 ) then
open(12, file='Simulated_Data/Simulated_data031.txt')
open(13, file='Simulated_Data/Simulated_HC011.txt')
open(15, file='Simulated_Data/Random_jzt011.txt')
open(16, file='Simulated_Data/Simulated_offers011.txt')

else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.32 ) then
open(12, file='Simulated_Data/Simulated_data032.txt')
open(13, file='Simulated_Data/Simulated_HC012.txt')
open(15, file='Simulated_Data/Random_jzt012.txt')
open(16, file='Simulated_Data/Simulated_offers012.txt')
end if 


do jfix=1, min(n_fix, I_state_fix) ; do je = 1 , n_e ; do jht=1, min(n_ht, I_state_ht)
    call TYPE1_DET(je,jht,jfix,jtype1)
do jsample=1,int(n_samp(jtype1))   
    
    ! read the same random draws as benchmark
      if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.GE.1) then
        do jno=1, 443
            read (14,*) ran_hypo(jno) ! 
        end do
      end if 
      
    ! read the HC and zt from saved files (experiments 9-12) if experiments are 17 to 20    
    if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.17) then
        do jno=1, 34
            read (13,*) HC_hypo_9(jno)
            read (15,*) zt_hypo_9(jno)
        end do
    else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.18) then
                do jno=1, 24
            read (13,*) HC_hypo_10(jno)
            read (15,*) zt_hypo_10(jno)
        end do
    else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.19) then
                        do jno=1, 14
            read (13,*) HC_hypo_11(jno)
            read (15,*) zt_hypo_11(jno)
        end do
    else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.20) then
                                do jno=1, 4
            read (13,*) HC_hypo_12(jno)
            read (15,*) zt_hypo_12(jno)
                                end do
   else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.27) then
                do jno=1, 24
            read (16,*) o_hypo_10(jno)
            end do 
            
 else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.29) then 
                            do jno=1, 34
            read (13,*) HC_hypo_9(jno)
            read (15,*) zt_hypo_9(jno)
            read (16,*) o_hypo_9(jno)
        end do
              
            
 else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.30) then
                     do jno=1, 24
            read (13,*) HC_hypo_10(jno)
            read (15,*) zt_hypo_10(jno)
            read (16,*) o_hypo_10(jno)
                     end do
                     
 else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.31) then 
                     do jno=1, 14
            read (13,*) HC_hypo_11(jno)
            read (15,*) zt_hypo_11(jno)
            read (16,*) o_hypo_11(jno)
                     end do
                     
else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.32) then 
                     do jno=1, 4
            read (13,*) HC_hypo_12(jno)
            read (15,*) zt_hypo_12(jno)
            read (16,*) o_hypo_12(jno)
        end do
        
    end if 
    
    
do jage=1,n_age 


! *************   FIRST PERIOD WHEN AGE=1    **************
if (jage.eq.1) then 
               jassets=assets_initial(je)
               jx=0 ! 0 years of experience
               jx_int=0.d0    
                    rb=ran0(idum) 
                      if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.0) then
                        write(14,1005) rb 
                        else  
                             rb=ran_hypo(1) 
                         end if 
                    bb=1
                    prob=In_h(je,bb,jht)
                 57 if (rb.lt.prob) then
                    jh=bb
                    else 
                    bb=bb+1
                    if (bb.le.n_h) then
		            prob=prob+In_h(je,bb,jht)
		            go to 57
		            else 
		            jh=n_h 
                    end if 
                    end if    
                    
                      jh_alt=jh
                    
!Pick jr from the initial distribution of health at age 25 from data
                    rb=ran0(idum) 
                      if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.0) then
                          write(14,1005) rb
                          else
                           rb=ran_hypo(2) 
                       end if 
                    bb=1
                    prob=In_r(je,bb)
                 58 if (rb.lt.prob) then
                    jr=bb
                    else 
                    bb=bb+1
                       if (bb.le.n_r) then
		            prob=prob+In_r(je,bb)
		            go to 58
                    else 
		            jr=n_r 
                    end if 
                    end if 
                    
                    jr_alt=jr 
                    
                       if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.6) then 
                       jr=1
                       end if
                    
!Pick jmar from the initial distribution at age 25 from data
                    rb=ran0(idum)   
                                          if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.0) then
                          write(14,1005) rb 
                          else 
                           rb=ran_hypo(3) 
                       end if 
                    prob=In_mar(je,jh)
                    if (rb.le.prob) then
                    jmar=2
                    else 
                    jmar=1
                    end if
                    
!draw spouse work status
                        rb=ran0(idum) 
                                          if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.0) then
                          write(14,1005) rb 
                          else
                           rb=ran_hypo(4) 
                       end if 
                    
                    if (jmar.eq.1) then
                        jos=1 ! assign any number, does not matter
                    else if (jmar.eq.2) then ! draw spouse work status
                        if (rb.lt.(prob_os(jage,je,jh,1,jfix))) then
                        jos=1
                        else 
                        jos=2
                        end if 
                    end if 
                    
!Draw employment offer
                  rb=ran0(idum) 
                                               if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.0) then
                          write(14,1005) rb 
                          else
                           rb=ran_hypo(5) 
                       end if 
                                  
		            bb=1
		            prob=Pr_o_25(bb, je)
85		            if (rb.lt.prob) then
		            jo=bb
		            else 
		            bb=bb+1
		            if (bb.le.n_o) then
		            prob=prob+Pr_o_25(bb, je)
		            go to 85
                    else 
                        jo=n_o
		            end if
                    end if
		            
    
                    rb=ran0(idum) ! draw a random # between 0-1
                    
                        if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.0) then
                          write(14,1005) rb 
                          else
                           rb=ran_hypo(6) 
                          end if 
                     
                         
		            bb=1
		            prob=p_zt(je,bb,2)
86		            if (rb.lt.prob) then
		            jzt=bb
		            else 
		            bb=bb+1
		            if (bb.le.n_zt) then
		            prob=prob+p_zt(je,bb,2)
		            go to 86
                    else 
                        jzt=n_zt
		            end if
                    end if 
                    

                                   
                    
end if  !for jage=1         
!*********************************************************

! *************   ALL AGES<N_RET    **************
if (jage.lt.n_ret) then	
    
        !experiments 
        if (jExperimentNumber.EQ.0.AND.(jBenchmark_exp.EQ.9).AND.jage.GE.7.AND.jage.LE.40 ) then 
        write(13,1004) jx_int
        write(15,1006) jzt
         write(16,1006) jo
    else if (jExperimentNumber.EQ.0.AND.(jBenchmark_exp.EQ.10).AND.jage.GE.17.AND.jage.LE.40 ) then 
        write(13,1004) jx_int
        write(15,1006) jzt
        write(16,1006) jo
    else if (jExperimentNumber.EQ.0.AND.(jBenchmark_exp.EQ.11).AND.jage.GE.27.AND.jage.LE.40 ) then 
        write(13,1004) jx_int
        write(15,1006) jzt
         write(16,1006) jo
    else if (jExperimentNumber.EQ.0.AND.(jBenchmark_exp.EQ.12).AND.jage.GE.37.AND.jage.LE.40 ) then 
        write(13,1004) jx_int
        write(15,1006) jzt
         write(16,1006) jo
    end if 
    
    !experiments 
        if (jExperimentNumber.EQ.0.AND.(jBenchmark_exp.EQ.17).AND.jage.GE.7.AND.jage.LE.40 ) then ! age 31 OR OLDER
            if (HC_hypo_9(jage-6).NE.0.0d0) then
        jx_int=HC_hypo_9(jage-6) 
            end if 
            if (zt_hypo_9(jage-6).NE.0) then
         jzt= zt_hypo_9(jage-6)
         end if 
        else if (jExperimentNumber.EQ.0.AND.(jBenchmark_exp.EQ.18).AND.jage.GE.17.AND.jage.LE.40 ) then 
                        if (HC_hypo_10(jage-16).NE.0.0d0) then
        jx_int=HC_hypo_10(jage-16) 
                        end if 
                        if (zt_hypo_10(jage-16).NE.0) then
          jzt= zt_hypo_10(jage-16)
                        end if 
                        
        else if (jExperimentNumber.EQ.0.AND.(jBenchmark_exp.EQ.27).AND.jage.GE.17.AND.jage.LE.40 ) then
                             if (o_hypo_10(jage-16).NE.0) then
        jo=o_hypo_10(jage-16) 
                        end if 
        else if (jExperimentNumber.EQ.0.AND.(jBenchmark_exp.EQ.19).AND.jage.GE.27.AND.jage.LE.40 ) then 
            if (HC_hypo_11(jage-26).NE.0.0d0) then
        jx_int=HC_hypo_11(jage-26)
            end if 
            if (zt_hypo_11(jage-26).NE.0) then
          jzt= zt_hypo_11(jage-26)
          end if 
        else if (jExperimentNumber.EQ.0.AND.(jBenchmark_exp.EQ.20).AND.jage.GE.37.AND.jage.LE.40 ) then 
            if (HC_hypo_12(jage-36).NE.0.0d0) then
         jx_int=HC_hypo_12(jage-36) 
            end if 
            if (zt_hypo_12(jage-36).NE.0) then
           jzt= zt_hypo_12(jage-36)
            end if
            
            
else if (jExperimentNumber.EQ.0.AND.(jBenchmark_exp.EQ.29).AND.jage.GE.7.AND.jage.LE.40 ) then ! age 31 OR OLDER
            if (HC_hypo_9(jage-6).NE.0.0d0) then
        jx_int=HC_hypo_9(jage-6) 
            end if 
            if (zt_hypo_9(jage-6).NE.0) then
         jzt= zt_hypo_9(jage-6)
            end if 
            if (o_hypo_9(jage-6).NE.0) then
        jo=o_hypo_9(jage-6) 
                        end if 
            
else if (jExperimentNumber.EQ.0.AND.(jBenchmark_exp.EQ.30).AND.jage.GE.17.AND.jage.LE.40 ) then
                                               if (HC_hypo_10(jage-16).NE.0.0d0) then
        jx_int=HC_hypo_10(jage-16) 
                        end if 
                        if (zt_hypo_10(jage-16).NE.0) then
          jzt= zt_hypo_10(jage-16)
                        end if
                             if (o_hypo_10(jage-16).NE.0) then
        jo=o_hypo_10(jage-16) 
                             end if 
                             
else if (jExperimentNumber.EQ.0.AND.(jBenchmark_exp.EQ.31).AND.jage.GE.27.AND.jage.LE.40 ) then                             
                        if (HC_hypo_11(jage-26).NE.0.0d0) then
        jx_int=HC_hypo_11(jage-26)
            end if 
            if (zt_hypo_11(jage-26).NE.0) then
          jzt= zt_hypo_11(jage-26)
            end if 
                                    if (o_hypo_11(jage-26).NE.0) then
        jo=o_hypo_11(jage-26) 
                             end if 
            
else if (jExperimentNumber.EQ.0.AND.(jBenchmark_exp.EQ.32).AND.jage.GE.37.AND.jage.LE.40 ) then 
               if (HC_hypo_12(jage-36).NE.0.0d0) then
         jx_int=HC_hypo_12(jage-36) 
            end if 
            if (zt_hypo_12(jage-36).NE.0) then
           jzt= zt_hypo_12(jage-36)
            end if
                                  if (o_hypo_12(jage-36).NE.0) then
        jo=o_hypo_12(jage-36) 
                             end if
                             
    end if 
    
    !Draw mortality shock            
                     rb=ran0(idum)
                                 if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.0) then
                                write(14,1005) rb 
                                end if 
                                 if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.GE.1) then
                                rb=ran_hypo(7+11*(jage-1))
                                    end if
		             if (rb.gt.p_surv(jh,je,jmar,  jage)) then !individual dies	 
     !randwom draws to files                     
if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.0) then
    if (jage.eq.40) then
        do jno=1,7
            rb=ran0(idum)
            write(14,1005) rb
        end do
        
        else 
                 do jno=1, 18 +(39-jage)*11 ! this makes sure that all individuals have 443 observations written to the random draws file- makes it easier to read.  
                     rb=ran0(idum)
            write(14,1005) rb
                 end do 
        end if 
end if      

       ! HC and zt shocks to be as in experiments 9-12
    !experiments 
        if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.9.AND.jage.le.39) then ! age 31 OR OLDER
            do jno=1,(40-max(jage,6))
       write(13,1004)  0.0d0
        write(15,1006)  0
        write(16,1006)  0
        end do
        else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.10.AND.jage.le.39) then 
            do jno=1,(40-max(jage,16))
       write(13,1004)  0.0d0
        write(15,1006)  0
         write(16,1006)  0
        end do
        else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.11.AND.jage.le.39) then 
            do jno=1,(40-max(jage,26))
       write(13,1004)  0.0d0
        write(15,1006)  0
        write(16,1006)  0
        end do
        else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.12.AND.jage.le.39) then 
             do jno=1,(40-max(jage,36))
            write(13,1004)  0.0d0
        write(15,1006)  0
        write(16,1006)  0
        end do
    end if 
    
    

                     go to 33 !simulate new individual 
                     end if

call TYPE4_DET(jmar,jos,jtype4) ! if jage is not 1, then jmar and jos have been already drawn towards the end of the code.                                            
call TYPE3_DET(jh,jr,jtype3) ! we need this index only for non-retirees

! determine optimal employment status "jo_act"
   !The individual decides whether to accept the offer or decline
   !interpolate with respect to jas and jx.

 if (jo.eq.1) then 
         jo_act=1  ! THIS IS YOUR ACTUAL EMPLOYMENT STATUS 
         
 else  ! HAVE TO DECIDE BETWEEN WORKING AND NOT WORKING
         EXPECTED_B(jo)=0.d0
         EXPECTED_B(1)=0.0d0

call interpol2d(as(1:n_as,jtype1,jage), xfn(   1: n_hc_grid , jage)  ,   B(1,jtype4,jzt,1:n_as,jtype3,jtype1,1: n_hc_grid,jage) , jassets, jx_int, EXPECTED_B(1))    
call interpol2d(as(1:n_as,jtype1,jage), xfn(   1: n_hc_grid , jage)  ,  B(jo,jtype4,jzt,1:n_as,jtype3,jtype1,1: n_hc_grid,jage) , jassets, jx_int, EXPECTED_B(jo))              
         If   (EXPECTED_B(jo).gt.EXPECTED_B(1)) then  
           jo_act=jo   
         else 
           jo_act=1 
         end if 
                
 end if 
 
!Draw health shocks (TYPE2)
                  rb=ran0(idum) ! draw a random # between 0-1 
                  
            if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.0) then
            write(14,1005) rb
            else 
                rb=ran_hypo(8+11*(jage-1))
            end if
                 
		            bb=1
		            prob=shocks_risk(je,bb, jh, jr, jage) 
		         90 if (rb.lt.prob) then
		            jtype2=bb
		            else 
		            bb=bb+1
		            if (bb.le.n_type2) then
		            prob=prob+shocks_risk(je,bb, jh, jr, jage)
		            go to 90
		            else 
		            jtype2=n_type2 
		            end if
                    end if
                    
if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.1 ) then
 if (jtype2.eq.3.OR.jtype2.eq.4.OR.jtype2.eq.7.OR.jtype2.eq.8) then 
     jtype2=jtype2-2
end if
else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.2 ) then
 if (jtype2.eq.5.OR.jtype2.eq.6.OR.jtype2.eq.7.OR.jtype2.eq.8) then 
     jtype2=jtype2-4
end if	            
else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.3 ) then
 if (jtype2.eq.2.OR.jtype2.eq.4.OR.jtype2.eq.6.OR.jtype2.eq.8) then
     jtype2=jtype2-1
end if		         
else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.4 ) then
 if (jtype2.eq.3.OR.jtype2.eq.4)   then 
     jtype2=jtype2-2
 else if (jtype2.eq.5.OR.jtype2.eq.6) then 
      jtype2=jtype2-4
else if (jtype2.eq.7.OR.jtype2.eq.8) then 
       jtype2=jtype2-6
end if 
else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.5 ) then
 jtype2=1
else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.9.AND.jage.EQ.6 ) then ! age 30, do not allow for d shocks
     if (jtype2.eq.5.OR.jtype2.eq.6.OR.jtype2.eq.7.OR.jtype2.eq.8) then 
     jtype2=jtype2-4
     end if	 	   
else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.10.AND.jage.EQ.16 ) then ! age 40, do not allow for d shocks
         if (jtype2.eq.5.OR.jtype2.eq.6.OR.jtype2.eq.7.OR.jtype2.eq.8) then 
     jtype2=jtype2-4
     end if	 
else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.11.AND.jage.EQ.26 ) then ! age 50, do not allow for d shocks
         if (jtype2.eq.5.OR.jtype2.eq.6.OR.jtype2.eq.7.OR.jtype2.eq.8) then 
     jtype2=jtype2-4
     end if	 
else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.12.AND.jage.EQ.36 ) then ! age 60, do not allow for d shocks
         if (jtype2.eq.5.OR.jtype2.eq.6.OR.jtype2.eq.7.OR.jtype2.eq.8) then 
     jtype2=jtype2-4
     end if	   
else if (jExperimentNumber.EQ.0.AND.(jBenchmark_exp.EQ.13.OR.jBenchmark_exp.EQ.17.OR.jBenchmark_exp.EQ.29).AND.jage.EQ.6 ) then ! age 30, ALL GET du shocks
     if (jtype2.eq.1.OR.jtype2.eq.2.OR.jtype2.eq.3.OR.jtype2.eq.4) then 
     jtype2=jtype2+4
     end if	   
else if (jExperimentNumber.EQ.0.AND.(jBenchmark_exp.EQ.14.OR.jBenchmark_exp.EQ.18.OR.jBenchmark_exp.EQ.27.OR.jBenchmark_exp.EQ.30).AND.jage.EQ.16 ) then ! age 40, ALL GET du shocks
      if (jtype2.eq.1.OR.jtype2.eq.2.OR.jtype2.eq.3.OR.jtype2.eq.4) then 
     jtype2=jtype2+4
     end if	 
else if (jExperimentNumber.EQ.0.AND.(jBenchmark_exp.EQ.15.OR.jBenchmark_exp.EQ.19.OR.jBenchmark_exp.EQ.31).AND.jage.EQ.26 ) then ! age 50, ALL GET du shocks
     if (jtype2.eq.1.OR.jtype2.eq.2.OR.jtype2.eq.3.OR.jtype2.eq.4) then 
     jtype2=jtype2+4
     end if	 
else if (jExperimentNumber.EQ.0.AND.(jBenchmark_exp.EQ.16.OR.jBenchmark_exp.EQ.20.OR.jBenchmark_exp.EQ.32).AND.jage.EQ.36 ) then ! age 60, ALL GET du shocks
      if (jtype2.eq.1.OR.jtype2.eq.2.OR.jtype2.eq.3.OR.jtype2.eq.4) then 
     jtype2=jtype2+4
     end if	 
end if	                    
                    
if (jtype2.eq.1) then
jdu=1 ; jdp=1 ; js=1
else if  (jtype2.eq.2) then
jdu=1 ; jdp=2 ; js=1
else if  (jtype2.eq.3) then
jdu=1 ; jdp=1 ; js=2
else if  (jtype2.eq.4) then
jdu=1 ; jdp=2 ; js=2
else if  (jtype2.eq.5) then
jdu=2 ; jdp=1 ; js=1
else if  (jtype2.eq.6) then
jdu=2 ; jdp=2 ; js=1
else if  (jtype2.eq.7) then
jdu=2 ; jdp=1 ; js=2
else if  (jtype2.eq.8) then
jdu=2 ; jdp=2 ; js=2
end if
		             
!draw the medical expenditure catastrophic shock
                    rb=ran0(idum) ! draw a random # between 0-1
                    
                                if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.0) then
            write(14,1005) rb 
            else

        rb=ran_hypo(9+11*(jage-1))
                      end if
                                     
		            bb=1
		            prob= Prob_MEcat(bb)
		         28 if (rb.lt.prob) then
		            jMEcat=bb
		            else 
		            bb=bb+1
		            if (bb.le.n_MEcat) then
		            prob=prob+Prob_MEcat(bb)
		            go to 28
		            else 
		            jMEcat= n_MEcat  
		            end if
                    end if  
                    
                    
                    ! draw whether you have the option to not pay and treat. 
                                        rb=ran0(idum) 
                                                                        if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.0) then
            write(14,1005) rb 
            else

        rb=ran_hypo(10+11*(jage-1))
                      end if
                                        
                    if ((jo_act.eq.3).OR.(jo_act.eq.5)) then
                        joption_pay=2 ! for those with ESHI, they all have option not to pay and get treated
                    else 
                    
            		prob= Shock_type(1,jo_act,jh,jage) + Shock_type(2,jo_act,jh,jage)   !Prob_option_pay(bb)
                  if (rb.le.prob) then
		            joption_pay=1 ! cannot treat and not pay
		            else 
		            joption_pay=2 ! can treat and not pay
                    end if  
                    end if 
                    
                    
                                        rb=ran0(idum)
                                          if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.0) then
                                            write(14,1005) rb 
                                            else
                         rb=ran_hypo(11+11*(jage-1))
                      end if
                      
                    if ((joption_pay.eq.1).and.(jtype2.ne.1)) then
                    if (rb.lt.(Shock_type(1,jo_act,jh,jage)/(Shock_type(1,jo_act,jh,jage)+Shock_type(2,jo_act,jh,jage)))) then ! you cannot treat 
                        juntreated=1 ! cannot treat
                    else 
                        juntreated=2 ! this means you have option to treat
                    end if 
                    else 
                        juntreated=2 ! set to 2 for everyone else 
                    end if 
                    
         if (oop(jtype2,jh,jage,jMEcat,jo_act).le.MC_TS(1)) then 
        jmc=1
        else if (oop(jtype2,jh,jage,jMEcat,jo_act).le.MC_TS(2)) then
        jmc=2
        else if (oop(jtype2,jh,jage,jMEcat,jo_act).le.MC_TS(3)) then
        jmc=3
        else
        jmc=4
        end if
 
!given the optimal work decision and sick days, determine human capital at the beginning of next period            

    if (jo_act.eq.1) then 
     jx_int_p=jx_int
     else if (jo_act.eq.2.OR.jo_act.eq.3) then  
     jx_int_p= min((jx_int+ (max( (hrs_pt - phi(je,jh,jtype2) ), 0.0d0) /hrs_ft )),  dble(n_hc) )
     else if (jo_act.eq.4.OR.jo_act.eq.5) then 
     jx_int_p=   min((jx_int+ (max((hrs_ft - phi(je,jh,jtype2) ), 0.0d0) /hrs_ft )),  dble(n_hc) ) ! jx_int goes up by 1 if you work full time and have no sick days. 
     end if
     
if (jage.le.n_ret-2) then    
do jzhc_p=1,n_zhc ! no longer a feature. no shocks. n_zhc=1
      HK_INTERPOL_sim(jzhc_p) =max(0.0d0,  min (jx_int_p *zhc(jzhc_p), dble(n_hc))) 
end do       
end if

!Find wages associated with current human capital level "jx_int"-- NEED TO INTERPOLATE: get "Z_INT" and "Z_INT_OFFER"
Z_INT_OFFER=0.0d0
Z_INT=0.0d0

call InterpolateScalar(HK(0 : n_hc), z(0 : n_hc,jh,je,jfix,jzt), jx_int, Z_INT_OFFER)
    if (jo_act.ge.2) then  
     Z_INT=Z_INT_OFFER
    end if
    
    
if (jage.eq.(n_ret-1)) then 
!Determine human capital group for the purpose of SS next period for those aged 64 (jhc_ret)
if (jx_int.le.hc_TS(1,je)) then 
    jhc_ret=1
else if (jx_int.le.hc_TS(2,je)) then 
    jhc_ret=2
else 
    jhc_ret=3
end if
else 
   jhc_ret=1  ! this will not matter for anything. 
end if 
 

! calculate incomes if pay or not pay
 actual_earn=0.0d0
 income_bt=0.0d0
 tax_SS=0.0d0
 tax_med=0.0d0
 
  actual_earn = max(0.0d0 , ( wage_penalty(jo_act,je)*Z_INT * (hrs_offer(jo_act) -phi(je,jh,jtype2) )) ) ! earnings net of sick days
  income_bt= ((interest-1.0d0))*jassets +  actual_earn  ! income before tax is equal to interest income plus earnings net of sick days
 
if ((jo_act.eq.3).OR.(jo_act.eq.5)) then ! husband holds insurance 
tax_SS = tau_ss*max(0.0d0, min(yhat ,  actual_earn -prem_ephi(jo_act,jtype4) ))   + tau_ss*max(0.0d0, min(yhat ,  inc_spouse(jage,je,jh,jtype4,jfix)  )) 
tax_med = tau_med* max(0.0d0,  ( actual_earn -prem_ephi(jo_act,jtype4))) +  tau_med* max(0.0d0,  ( inc_spouse(jage,je,jh,jtype4,jfix) ))
else ! subtract from wife. 
tax_SS = tau_ss*max(0.0d0, min(yhat ,  actual_earn  ))   + tau_ss*max(0.0d0, min(yhat ,  inc_spouse(jage,je,jh,jtype4,jfix) -prem_ephi(jo_act,jtype4) )) 
tax_med = tau_med* max(0.0d0,  ( actual_earn )) +  tau_med* max(0.0d0,  ( inc_spouse(jage,je,jh,jtype4,jfix) -prem_ephi(jo_act,jtype4)))
end if   
  
 ! for all of these, pay OOP =1, not pay OOP=2 
 income_pay(:)=0.0d0
 Total_tax_pay(:)=0.0d0
 base_2_pay(:)=0.0d0
 base_1_pay(:)=0.0d0
 
        ! if pay OOP 
  income_pay(1) = max(0.d0,  income_bt + inc_spouse(jage,je,jh,jtype4,jfix) -prem_ephi(jo_act,jtype4) -tax_SS - tax_med - max(0.d0, oop(jtype2,jh,jage,jMEcat,jo_act) + oop_spouse(jage,je,jtype4,jo_act) - 0.075d0*(income_bt+  inc_spouse(jage,je,jh,jtype4,jfix)) ) )  ! taxable income
  Total_tax_pay(1) = max(0.0d0,  (( (income_pay(1)*5200.0d0) - tax_lambda(jmar) *( (income_pay(1)*5200.0d0)**(1.0d0-tax_tau(jmar))) )/5200 )) +  tax_SS + tax_med  ! tax you have to pay
  base_2_pay(1) =  prem_ephi(jo_act,jtype4) +  oop(jtype2,jh,jage,jMEcat,jo_act) +  oop_spouse(jage,je,jtype4,jo_act) +Total_tax_pay(1)  ! total resources subtracted 
  base_1_pay(1) =  interest*jassets + inc_spouse(jage,je,jh,jtype4,jfix)  + actual_earn - base_2_pay(1)        ! cash on hand (if you pay the OOP)
      ! if not pay OOP
  income_pay(2) = max(0.d0,  income_bt + inc_spouse(jage,je,jh,jtype4,jfix) -prem_ephi(jo_act,jtype4)  -tax_SS - tax_med -  max(0.d0,  oop_spouse(jage,je,jtype4,jo_act) - 0.075d0*(income_bt+  inc_spouse(jage,je,jh,jtype4,jfix)) ) )  ! taxable income
  Total_tax_pay(2) = max(0.0d0,  (( (income_pay(2)*5200.0d0) - tax_lambda(jmar) *( (income_pay(2)*5200.0d0)**(1.0d0-tax_tau(jmar))) )/5200 )) +  tax_SS + tax_med  ! tax you have to pay
  base_2_pay(2) =  prem_ephi(jo_act,jtype4) +  oop_spouse(jage,je,jtype4,jo_act) + Total_tax_pay(2)  ! total resources subtracted 
  base_1_pay(2) =  interest*jassets + inc_spouse(jage,je,jh,jtype4,jfix)  + actual_earn - base_2_pay(2)        ! cash on hand (if you do not pay the OOP)
  
  
  
!Determine income group
          if (income_bt.le.INC_TS(1)) then 
        jinc=1
        else if (income_bt.le.INC_TS(2)) then 
        jinc=2
        else if (income_bt.le.INC_TS(3)) then 
        jinc=3
        else if (income_bt.le.INC_TS(4)) then 
        jinc=4
        else 
        jinc=5
        end if
    

        
!DETERMINE IF THE PERSON WILL PAY/TREAT and OPTIMAL CONSUMPTION/SAV
!(determine assets next period): jcons_act, jas_p_act , jsav_act

  !initiate        
CON_sim(:)=0.0d0
  ASP_sim(:)=0.d0
  B_sim(:)=0.d0
  TR_sim(:)=0.0d0
  
  CON=0.0d0
  ASP=0.d0
  UTI=0.d0
  TR=0.0d0
  MAX_SAV=0.0d0 
  V_cf_trpay=0.0d0
  V_c_all_trpay=0.0d0
  V_death=0.0d0
        jcons_act=0.0d0
        jas_p_act=0.0d0
        jtr_act=0.0d0
        B_act=0.0d0
  
  
!If I don't have a shock, I pay medical bill. 
if (jtype2.eq.1) then
     jtrpay=1
    I_treat=2
    I_pay=1
    
    !If on cons floor
    if ( base_1_pay(1) .le. (cons_floor_at(jage, jmar,je,jh))) then 
    jcons_act=cons_floor(jage, jmar,je,jh)  ! cons floor is at the HH level; CON is after tax HH consumption - it's what goes into utility fn
    jas_p_act=0.d0 
    jtr_act=  cons_floor_at(jage, jmar,je,jh) - base_1_pay(1)
    
    ! still need to calculate V because we need it all the time for welfare analysis
    ! get utility
    call func_uti(jage,1,je,jh,jo,jtype4,jtype2,jcons_act,UTI) ! second argument is trpay which we set to 1
             V_death =  Bequest_val(0.0d0)  + cost_death
         !Calculate expected value associated with MAX_SAV assets tomorrow
         EXPECTED_B_PRIME=0.0d0
         
         if (jage.eq.(n_ret-1)) then 

                    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1, n_MEcat; do jmar_p=1, n_mar
                     EXPECTED_B_PRIME = EXPECTED_B_PRIME + Tp_h(jh_p, jh,je,jht,jage, jtype2,2,1)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p)* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (RET_PRIME(jhc_ret,jMEcat_p,1,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p) * p_surv(jh_p,je,jmar_p, jage+1) +  V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)) )                                      
                    end do; end do; end do  ; end do; end do
  
         else 
             EXP_L_prime_INT(:,:, :,:,:)=0.0d0
                    do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar;  do jzhc_p=1,n_zhc
                    call TYPE3_DET(jh_p,jr_p,jtype3_p) 
      !interpolate with respect to human capital. Assets will be 0 next period.           
     if (jo_act.eq.1) then
     call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,1), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )
     else
     call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,2), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )
     end if 
     
     EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(1,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))
                     end do ; end do ; end do   ; end do; end do  
  
          end if 
                 
         B_sim(1) = UTI + beta_hat(je)* EXPECTED_B_PRIME 

         
    
    else
              jtr_act=0.d0 
! find maximum savings on the grid (jsav_upper)
  MAX_SAV= income_bt + inc_spouse(jage,je,jh,jtype4,jfix) -base_2_pay(1)  - cons_floor_at(jage, jmar,je,jh)   
  
  !find lower bound 
          jsav_upper=0
          do while ( MAX_SAV .ge. sav_sim(jsav_upper+1) )
           jsav_upper=jsav_upper +1
           if (jsav_upper.EQ.n_sav_sim) exit
          end do 
          
               if (jsav_upper.lt.1) then
                 PRINT *,'Error'
               end if
               
!VALUE ASSOCIATED WITH CONSUMING on the consumption floor and saving max amount: V_cf_trpay
         call func_uti(jage,1,je,jh,jo,jtype4,jtype2,cons_floor(jage, jmar,je,jh),UTI) 
         V_death =  Bequest_val(MAX_SAV+jassets)  + cost_death
         !Calculate expected value associated with MAX_SAV assets tomorrow
         EXPECTED_B_PRIME=0.0d0
         
         if (jage.eq.(n_ret-1)) then 

                    RET_PRIME_INT(:,:,:,:,:)=0.0d0
                    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1, n_MEcat; do jmar_p=1, n_mar
                    !interpolate for assets
                        call InterpolateScalar(as(1:n_as,jtype1,jage+1), RET_PRIME(jhc_ret,jMEcat_p,1:n_as,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p), (MAX_SAV+jassets), RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) )
                     EXPECTED_B_PRIME = EXPECTED_B_PRIME +Tp_h(jh_p, jh,je,jht,jage, jtype2,2,1)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p)* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p)* p_surv(jh_p,je,jmar_p, jage+1) +  V_death * (1.0-p_surv(jh_p, je,jmar_p, jage+1)) )  
                    end do; end do; end do  ; end do; end do
  
         else 
             EXP_L_prime_INT(:,:, :,:,:)=0.0d0
                    do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar;  do jzhc_p=1,n_zhc
                    call TYPE3_DET(jh_p,jr_p,jtype3_p) 
      !interpolate with respect to human capital and assets   
                    ! the value function tomorrow depends on whether I'm working today
                    if (jo_act.eq.1) then
                    call interpol2d(as(1:n_as,jtype1,jage+1), xfn(   1: n_hc_grid , jage+1), EXP_L_prime(1:n_as,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,1), (MAX_SAV+jassets), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p))                           
                    else 
                    call interpol2d(as(1:n_as,jtype1,jage+1), xfn(   1: n_hc_grid , jage+1), EXP_L_prime(1:n_as,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,2), (MAX_SAV+jassets), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p))                           
                    end if 
                    
                    EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(1,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))
                     end do ; end do ; end do   ; end do; end do  
  
         end if 
                 
         V_cf_trpay = UTI + beta_hat(je)* EXPECTED_B_PRIME  

         
 !VALUE ASSOCIATED WITH CONSUMING ALL AND 0 ASSETS: V_c_all_trpay
         
           call func_uti(jage,1,je,jh,jo,jtype4,jtype2,(base_1_pay(1)/(1.0d0+taxc)),UTI) 
         V_death =  Bequest_val(0.0d0)  + cost_death
         !Calculate expected value associated with MAX_SAV assets tomorrow
         EXPECTED_B_PRIME=0.0d0
         
         if (jage.eq.(n_ret-1)) then 

                    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1, n_MEcat; do jmar_p=1, n_mar
                     EXPECTED_B_PRIME = EXPECTED_B_PRIME + Tp_h(jh_p, jh,je,jht,jage, jtype2,2,1)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p)* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (RET_PRIME(jhc_ret,jMEcat_p,1,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p) * p_surv(jh_p,je,jmar_p, jage+1) +  V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)) )                                      
                    end do; end do; end do  ; end do; end do
  
         else 
             EXP_L_prime_INT(:,:, :,:,:)=0.0d0
                    do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar;  do jzhc_p=1,n_zhc
                    call TYPE3_DET(jh_p,jr_p,jtype3_p) 
      !interpolate with respect to human capital. Assets will be 0 next period.           
     if (jo_act.eq.1) then
     call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,1), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )
     else
     call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,2), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )
     end if 
     
     EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(1,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))
                     end do ; end do ; end do   ; end do; end do  
  
          end if 
                 
         V_c_all_trpay = UTI + beta_hat(je)* EXPECTED_B_PRIME 
         
          
!Iterate over possible savings and find optimal   
         maxB =V_cf_trpay
         jas_p_act=MAX_SAV+jassets
         jcons_act=cons_floor(jage, jmar,je,jh)
do jsav= jsav_upper , 1, -1
! get ASP, CON, and UTI associated with this jsav    
        ASP=  jassets+  sav_sim(jsav)
     if (ASP.le.0.d0)  exit   ! I already did the case where you consume all and save 0, so that's why we have "le," not "lt."
        V_death= cost_death + Bequest_val(ASP)  
        CON= ( income_bt + inc_spouse(jage,je,jh,jtype4,jfix) - base_2_pay(1) -sav_sim(jsav) )/(1.0d0+taxc)
        call  func_uti(jage,1,je,jh,jo_act,jtype4,jtype2,CON,UTI)           
 
!Solve continuation value given this jsav and jtrpay
          EXPECTED_B_PRIME=0.0d0
         if (jage.eq.(n_ret-1)) then !need to interpolate only for assets
                    RET_PRIME_INT(:,:,:,:,:)=0.0d0


                    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1, n_MEcat; do jmar_p=1, n_mar
                        call InterpolateScalar(as(1:n_as,jtype1,jage+1), RET_PRIME(jhc_ret,jMEcat_p,1:n_as,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p), ASP, RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) )
                     EXPECTED_B_PRIME = EXPECTED_B_PRIME +Tp_h(jh_p, jh,je,jht,jage, jtype2,2,1)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p)* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p)* p_surv(jh_p,je,jmar_p, jage+1) +  V_death * (1.0-p_surv(jh_p, je,jmar_p, jage+1)) )                                      
                    end do; end do; end do  ; end do; end do
  
         else 
             EXP_L_prime_INT(:,:, :,:,:)=0.0d0
                    do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar;  do jzhc_p=1,n_zhc
                    call TYPE3_DET(jh_p,jr_p,jtype3_p) 
                if (jo_act.eq.1) then
                    call interpol2d(as(1:n_as,jtype1,jage+1), xfn(   1: n_hc_grid , jage+1), EXP_L_prime(1:n_as,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,1), ASP, HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p))                           
               else 
                    call interpol2d(as(1:n_as,jtype1,jage+1), xfn(   1: n_hc_grid , jage+1), EXP_L_prime(1:n_as,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,2), ASP, HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p))                           
               end if 
               
                  EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(1,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act)  *(EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p)* p_surv(jh_p,je,jmar_p, jage+1) + V_death * (1.0-p_surv(jh_p,je,jmar_p, jage+1)) ) 
                    end do ; end do ; end do   ; end do; end do  
  
          end if 
                                 
B_act=UTI+ beta_hat(je)*EXPECTED_B_PRIME 

     if (B_act.lt.maxB) exit
              maxB=B_act  
              jas_p_act=ASP
              jcons_act=CON
    
       
51 end do !jsav
! Found jcons_act , jas_p_act and B_act assocaited with optimal savings. 
   
! now also need to compare with:    V_c_all_trpay. And if this is max, then set jas and jcons accordingly.
  if  (V_c_all_trpay.gt.maxB) then
      maxB =V_c_all_trpay
      jas_p_act=0.0d0
      jcons_act= base_1_pay(I_pay)/(1.0d0+taxc)
  end if 
        B_sim(1)=   maxB    
    end if 
       
    
 !If I have a shock   
else 
if (( base_1_pay(2) .le. (cons_floor_at(jage, jmar,je,jh))).AND.(juntreated.EQ.1).AND.(I_treat_Mcaid.EQ.0)) then ! if I'm on the cons floor even if I don't pay AND cannot treat
    

        jtrpay=3 ! you don't treat
        I_treat=1 ! not treat
        I_pay=2 ! not pay
    jcons_act=cons_floor(jage, jmar,je,jh)  ! cons floor is at the HH level; CON is after tax HH consumption - it's what goes into utility fn
    jas_p_act=0.d0 
    jtr_act=  cons_floor_at(jage, jmar,je,jh) - base_1_pay(2)
  

    call func_uti(jage,jtrpay,je,jh,jo,jtype4,jtype2,jcons_act,UTI) ! second argument is trpay which we set to 3
             V_death =  Bequest_val(0.0d0)  + cost_death
         !Calculate expected value associated with MAX_SAV assets tomorrow
         EXPECTED_B_PRIME=0.0d0
         
         if (jage.eq.(n_ret-1)) then 

                    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1, n_MEcat; do jmar_p=1, n_mar
                     EXPECTED_B_PRIME = EXPECTED_B_PRIME + Tp_h(jh_p, jh,je,jht,jage, jtype2,1,1)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p)* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (RET_PRIME(jhc_ret,jMEcat_p,1,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p) * p_surv(jh_p,je,jmar_p, jage+1) +  V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)) )                                      
                    end do; end do; end do  ; end do; end do
  
         else 
             EXP_L_prime_INT(:,:, :,:,:)=0.0d0
                    do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar;  do jzhc_p=1,n_zhc
                    call TYPE3_DET(jh_p,jr_p,jtype3_p) 
      !interpolate with respect to human capital. Assets will be 0 next period.           
     if (jo_act.eq.1) then
     call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,1), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )
     else
     call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,2), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )
     end if 
     
     EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(2,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))
                     end do ; end do ; end do   ; end do; end do  
  
          end if 
                 
         B_sim(jtrpay) = UTI + beta_hat(je)* EXPECTED_B_PRIME 
    
    
else if (( base_1_pay(2) .le. (cons_floor_at(jage, jmar,je,jh))).AND.((juntreated.EQ.2).OR.(I_treat_Mcaid.EQ.1))) then ! if I'm on the cons floor even if I don't pay AND can treat - then you always "treat and pay"
    jtrpay=1
    I_treat=2
    I_pay=1
    jcons_act=cons_floor(jage, jmar,je,jh)  ! cons floor is at the HH level; CON is after tax HH consumption - it's what goes into utility fn
    jas_p_act=0.d0 
    jtr_act=  cons_floor_at(jage, jmar,je,jh) - base_1_pay(1)
    
    ! get utility
    call func_uti(jage,jtrpay,je,jh,jo,jtype4,jtype2,jcons_act,UTI) ! second argument is trpay which we set to 1
             V_death =  Bequest_val(0.0d0)  + cost_death
         !Calculate expected value associated with MAX_SAV assets tomorrow
         EXPECTED_B_PRIME=0.0d0
         
         if (jage.eq.(n_ret-1)) then 

                    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1, n_MEcat; do jmar_p=1, n_mar
                     EXPECTED_B_PRIME = EXPECTED_B_PRIME + Tp_h(jh_p, jh,je,jht,jage, jtype2,2,1)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p)* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (RET_PRIME(jhc_ret,jMEcat_p,1,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p) * p_surv(jh_p,je,jmar_p, jage+1) +  V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)) )                                      
                    end do; end do; end do  ; end do; end do
  
         else 
             EXP_L_prime_INT(:,:, :,:,:)=0.0d0
                    do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar;  do jzhc_p=1,n_zhc
                    call TYPE3_DET(jh_p,jr_p,jtype3_p) 
      !interpolate with respect to human capital. Assets will be 0 next period.           
     if (jo_act.eq.1) then
     call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,1), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )
     else
     call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,2), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )
     end if 
     
     EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(1,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))
                     end do ; end do ; end do   ; end do; end do  
  
          end if 
                 
         B_sim(jtrpay) = UTI + beta_hat(je)* EXPECTED_B_PRIME 
    
else
    
if (joption_pay.eq.2) then   
do jtrpay=1, 3
    

    
!jtrpay=1 	pay and get treated
!jtrpay=2	don't pay, get treated
!jtrpay=3	don't pay, not treated
    
    if (jtrpay.eq.1) then
        I_treat=2 ! treat
        I_pay=1 ! pay
    else if (jtrpay.eq.2) then
        I_treat=2 ! treat
        I_pay=2 ! not pay
    else 
        I_treat=1 ! not treat
        I_pay=2 ! not pay
    end if 

!initiate
  CON=0.0d0
  ASP=0.d0
  UTI=0.d0
  TR=0.0d0
  MAX_SAV=0.0d0 
        jcons_act=0.0d0
        jas_p_act=0.0d0
        jtr_act=0.0d0
        B_act=0.0d0
        
!If you are on the consumption floor        
if ( base_1_pay(I_pay) .le. (cons_floor_at(jage, jmar,je,jh))) then ! you get transfer. you only get to this point if you are not on transfer if you don't pay. so you'll still have to calculate the value if you don't pay
 jcons_act=cons_floor(jage, jmar,je,jh)  ! cons floor is at the HH level; CON is after tax HH consumption - it's what goes into utility fn
 jas_p_act=0.d0 
 jtr_act=  cons_floor_at(jage, jmar,je,jh) - base_1_pay(I_pay) 
 call  func_uti(jage,jtrpay,je,jh,jo_act,jtype4,jtype2,cons_floor(jage, jmar,je,jh),UTI)
 V_death =  Bequest_val(0.0d0)  + cost_death 
 
!continuation value
 EXPECTED_B_PRIME=0.0d0
         if (jage.eq.(n_ret-1)) then 


                    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1, n_MEcat; do jmar_p=1, n_mar
                     EXPECTED_B_PRIME = EXPECTED_B_PRIME + Tp_h(jh_p, jh,je,jht,jage, jtype2,I_treat,jmc)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p)* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (RET_PRIME(jhc_ret,jMEcat_p,1,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p) * p_surv(jh_p,je,jmar_p, jage+1) +  V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)) )                                      
                    end do; end do; end do  ; end do; end do
  
         else 
             EXP_L_prime_INT(:,:, :,:,:)=0.0d0
                    do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar;  do jzhc_p=1,n_zhc
                    call TYPE3_DET(jh_p,jr_p,jtype3_p) 
      !interpolate with respect to human capital. Assets will be 0 next period.     
                    if (jo_act.eq.1) then 
     call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,1), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )
                    else 
     call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,2), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )

                    end if 
                    
     if ((jtrpay.eq.1).OR.(jtrpay.eq.2)) then ! treat
     EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(1,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))
      else               
          EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(jmc+1,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))
      end if 
      
     end do ; end do ; end do   ; end do; end do  
  
         end if 
         
         
         
                                 
maxB=UTI+ beta_hat(je)*EXPECTED_B_PRIME 
 
else
  jtr_act=0.d0
! find maximum savings on the grid (jsav_upper)
  MAX_SAV= income_bt + inc_spouse(jage,je,jh,jtype4,jfix) -base_2_pay(I_pay)  - cons_floor_at(jage, jmar,je,jh)   
  
  !find lower bound 
          jsav_upper=0
          do while ( MAX_SAV .ge. sav_sim(jsav_upper+1) )
           jsav_upper=jsav_upper +1
           if (jsav_upper.EQ.n_sav_sim) exit
          end do 
          
               if (jsav_upper.lt.1) then
                 PRINT *,'Error'
               end if 
                    
               
  !VALUE ASSOCIATED WITH CONSUMING on the consumption floor and saving max amount: V_cf_trpay
         call func_uti(jage,jtrpay,je,jh,jo,jtype4,jtype2,cons_floor(jage, jmar,je,jh),UTI) 
         V_death =  Bequest_val(MAX_SAV+jassets)  + cost_death
         !Calculate expected value associated with MAX_SAV assets tomorrow
         EXPECTED_B_PRIME=0.0d0
         
         if (jage.eq.(n_ret-1)) then 

                    RET_PRIME_INT(:,:,:,:,:)=0.0d0
                    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1, n_MEcat; do jmar_p=1, n_mar
                    !interpolate for assets
                        call InterpolateScalar(as(1:n_as,jtype1,jage+1), RET_PRIME(jhc_ret,jMEcat_p,1:n_as,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p), (MAX_SAV+jassets), RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) )
                     EXPECTED_B_PRIME = EXPECTED_B_PRIME +Tp_h(jh_p, jh,je,jht,jage, jtype2,I_treat,jmc)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p)* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p)* p_surv(jh_p,je,jmar_p, jage+1) + V_death * (1.0-p_surv(jh_p, je,jmar_p, jage+1)) )  
                    end do; end do; end do  ; end do; end do
  
         else 
             EXP_L_prime_INT(:,:, :,:,:)=0.0d0
                    do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar;  do jzhc_p=1,n_zhc
                    call TYPE3_DET(jh_p,jr_p,jtype3_p) 
      !interpolate with respect to human capital and assets  
                    if (jo_act.eq.1) then
                    call interpol2d(as(1:n_as,jtype1,jage+1), xfn(   1: n_hc_grid , jage+1), EXP_L_prime(1:n_as,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,1), (MAX_SAV+jassets), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p))                           
                    else 
                    call interpol2d(as(1:n_as,jtype1,jage+1), xfn(   1: n_hc_grid , jage+1), EXP_L_prime(1:n_as,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,2), (MAX_SAV+jassets), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p))                           

                    end if 
                    
      if ((jtrpay.eq.1).OR.(jtrpay.eq.2)) then ! treat              
                    EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(1,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))
      else
                              EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(jmc+1,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))

      end if 
      
      
      end do ; end do ; end do   ; end do; end do  
  
         end if 
                 
         V_cf_trpay = UTI + beta_hat(je)* EXPECTED_B_PRIME  

         
 !VALUE ASSOCIATED WITH CONSUMING ALL AND 0 ASSETS: V_c_all_trpay
         
           call func_uti(jage,jtrpay,je,jh,jo,jtype4,jtype2,(base_1_pay(I_pay)/(1.0d0+taxc)),UTI) 
         V_death =  Bequest_val(0.0d0)  + cost_death
         !Calculate expected value associated with MAX_SAV assets tomorrow
         EXPECTED_B_PRIME=0.0d0
         
         if (jage.eq.(n_ret-1)) then 
 
                    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1, n_MEcat; do jmar_p=1, n_mar
                     EXPECTED_B_PRIME = EXPECTED_B_PRIME + Tp_h(jh_p, jh,je,jht,jage, jtype2,I_treat,jmc)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p)* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (RET_PRIME(jhc_ret,jMEcat_p,1,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p) * p_surv(jh_p,je,jmar_p, jage+1) +  V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)) )                                      
                    end do; end do; end do  ; end do; end do
  
         else 
             EXP_L_prime_INT(:,:, :,:,:)=0.0d0
                    do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar;  do jzhc_p=1,n_zhc
                    call TYPE3_DET(jh_p,jr_p,jtype3_p) 
      !interpolate with respect to human capital. Assets will be 0 next period.   
                    if (jo_act.eq.1) then
     call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,1), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )
                    else 
     call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,2), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )

                    end if 
                    
     if ((jtrpay.eq.1).OR.(jtrpay.eq.2)) then ! treat 
     EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(1,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))
     
     else
              EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(jmc+1,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))

     end if 
     
     
     end do ; end do ; end do   ; end do; end do  
  
          end if 
                 
         V_c_all_trpay = UTI + beta_hat(je)* EXPECTED_B_PRIME 
         
          
!Iterate over possible savings and find optimal   
         maxB =V_cf_trpay
         jas_p_act=MAX_SAV+jassets
         jcons_act=cons_floor(jage, jmar,je,jh)
do jsav= jsav_upper , 1, -1
! get ASP, CON, and UTI associated with this jsav    
        ASP=  jassets+  sav_sim(jsav)
     if (ASP.le.0.d0)  exit   ! I already did the case where you consume all and save 0, so that's why we have "le," not "lt."
        V_death= cost_death + Bequest_val(ASP)  
        CON= ( income_bt + inc_spouse(jage,je,jh,jtype4,jfix) - base_2_pay(I_pay) -sav_sim(jsav) )/(1.0d0+taxc)
        call  func_uti(jage,jtrpay,je,jh,jo_act,jtype4,jtype2,CON,UTI)           
 
!Solve continuation value given this jsav and jtrpay
          EXPECTED_B_PRIME=0.0d0
         if (jage.eq.(n_ret-1)) then !need to interpolate only for assets; for human capital, we discretize it into bigger categories, so no need to interpolate.
                    RET_PRIME_INT(:,:,:,:,:)=0.0d0


                    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1, n_MEcat; do jmar_p=1, n_mar
                        call InterpolateScalar(as(1:n_as,jtype1,jage+1), RET_PRIME(jhc_ret,jMEcat_p,1:n_as,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p), ASP, RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) )
                     EXPECTED_B_PRIME = EXPECTED_B_PRIME +Tp_h(jh_p, jh,je,jht,jage, jtype2,I_treat,jmc)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p)* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p)* p_surv(jh_p,je,jmar_p, jage+1) +  V_death * (1.0-p_surv(jh_p, je,jmar_p, jage+1)) )                                      
                    end do; end do; end do  ; end do; end do
  
         else 
             EXP_L_prime_INT(:,:, :,:,:)=0.0d0
                    do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar;  do jzhc_p=1,n_zhc
                    call TYPE3_DET(jh_p,jr_p,jtype3_p) 
                                        if (jo_act.eq.1) then

                    call interpol2d(as(1:n_as,jtype1,jage+1), xfn(   1: n_hc_grid , jage+1), EXP_L_prime(1:n_as,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,1), ASP, HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p))                           
                                        else 
                    call interpol2d(as(1:n_as,jtype1,jage+1), xfn(   1: n_hc_grid , jage+1), EXP_L_prime(1:n_as,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,2), ASP, HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p))                           

                                        end if 
                                        
                                     
                    
               if ((jtrpay.eq.1).OR.(jtrpay.eq.2)) then ! treat     
                    EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(1,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act)  *(EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p)* p_surv(jh_p,je,jmar_p, jage+1) + V_death * (1.0-p_surv(jh_p,je,jmar_p, jage+1)) ) 
               else 
                                       EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(jmc+1,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act)  *(EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p)* p_surv(jh_p,je,jmar_p, jage+1) + V_death * (1.0-p_surv(jh_p,je,jmar_p, jage+1)) ) 

               end if 
               
               
               end do ; end do ; end do   ; end do; end do  
  
          end if 
                                 
B_act=UTI+ beta_hat(je)*EXPECTED_B_PRIME 


     if (B_act.lt.maxB) exit
              maxB=B_act  
              jas_p_act=ASP
              jcons_act=CON
    
       
 end do !jsav
! Found jcons_act , jas_p_act and B_act assocaited with optimal savings. 
   
! now also need to compare with:    V_c_all_trpay. And if this is max, then set jas and jcons accordingly.
  if  (V_c_all_trpay.gt.maxB) then
      maxB =V_c_all_trpay
      jas_p_act=0.0d0
      jcons_act= base_1_pay(I_pay)/(1.0d0+taxc)
  end if 
  

end if ! for being on the consumption floor

! save CONS, ASP, TR and B for the 3 jtrpay options. 

  CON_sim(jtrpay)=jcons_act
  ASP_sim(jtrpay)=jas_p_act
  TR_sim(jtrpay)=jtr_act
  B_sim(jtrpay)=maxB
  
end do ! jtrpay

! determine best option and reset I_treat and I_pay
if ((B_sim(1).eq.0.0d0).OR.(B_sim(2).eq.0.0d0).or.(B_sim(3).eq.0.0d0)) then
PRINT *,'Error'
pause
end if 
if ((B_sim(1).gt.B_sim(2)).AND.(B_sim(1).ge.B_sim(3))) then
    jtrpay = 1 
else if ((B_sim(2).ge.B_sim(1)).AND.(B_sim(2).ge.B_sim(3))) then
    jtrpay = 2 
else if ((B_sim(2).eq.B_sim(1)).AND.(B_sim(2).eq.B_sim(3))) then
    jtrpay = 2 
    else
    jtrpay = 3 
end if 

jtr_act= TR_sim(jtrpay)
jas_p_act=ASP_sim(jtrpay)
jcons_act=CON_sim(jtrpay)

 else if (((joption_pay.eq.1).AND.(juntreated.EQ.2)).OR.(((joption_pay.eq.1).AND.(juntreated.EQ.1)).AND.(I_treat_Mcaid.EQ.1).AND.(base_1_pay(1) .le. (cons_floor_at(jage, jmar,je,jh))) )) then ! if you don't have the option to treat and not pay. BUT HAVE OPTION TO TREAT

    do jtrpay=1, 3, 2 
    
! I_treat and I_pay    
    
!jtrpay=1 	pay and get treated
!jtrpay=2	don't pay, get treated
!jtrpay=3	don't pay, not treated
    
    if (jtrpay.eq.1) then
        I_treat=2 ! treat
        I_pay=1 ! pay
    else 
        I_treat=1 ! not treat
        I_pay=2 ! not pay
    end if 

!initiate
  CON=0.0d0
  ASP=0.d0
  UTI=0.d0
  TR=0.0d0
  MAX_SAV=0.0d0 
        jcons_act=0.0d0
        !jsav_act=0.0d0 
        jas_p_act=0.0d0
        jtr_act=0.0d0
        B_act=0.0d0
        
!If you are on the consumption floor        
if ( base_1_pay(I_pay) .le. (cons_floor_at(jage, jmar,je,jh))) then ! you get transfer. you only get to this point if you are not on transfer if you don't pay. so you'll still have to calculate the value if you don't pay
 jcons_act=cons_floor(jage, jmar,je,jh)  ! cons floor is at the HH level; CON is after tax HH consumption - it's what goes into utility fn
 jas_p_act=0.d0 
 jtr_act=  cons_floor_at(jage, jmar,je,jh) - base_1_pay(I_pay) 
 call  func_uti(jage,jtrpay,je,jh,jo_act,jtype4,jtype2,cons_floor(jage, jmar,je,jh),UTI)
 V_death =  Bequest_val(0.0d0)  + cost_death 
 
!continuation value
 EXPECTED_B_PRIME=0.0d0
         if (jage.eq.(n_ret-1)) then 
 

                    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1, n_MEcat; do jmar_p=1, n_mar
                     EXPECTED_B_PRIME = EXPECTED_B_PRIME + Tp_h(jh_p, jh,je,jht,jage, jtype2,I_treat,jmc)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p)* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (RET_PRIME(jhc_ret,jMEcat_p,1,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p) * p_surv(jh_p,je,jmar_p, jage+1) +  V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)) )                                      
                    end do; end do; end do  ; end do; end do
  
         else 
             EXP_L_prime_INT(:,:, :,:,:)=0.0d0
                    do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar;  do jzhc_p=1,n_zhc
                    call TYPE3_DET(jh_p,jr_p,jtype3_p) 
      !interpolate with respect to human capital. Assets will be 0 next period.  
                    if (jo_act.eq.1) then
     call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,1), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )
                    else 
     call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,2), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )

                    end if 
                    
     if ((jtrpay.eq.1).OR.(jtrpay.eq.2)) then ! treat 
     EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(1,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))
     else 
              EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(jmc+1,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))
     end if 
     
     
     end do ; end do ; end do   ; end do; end do  
  
          end if 
                                 
maxB=UTI+ beta_hat(je)*EXPECTED_B_PRIME 
 
else
  jtr_act=0.d0
! find maximum savings on the grid (jsav_upper)
  MAX_SAV= income_bt + inc_spouse(jage,je,jh,jtype4,jfix) -base_2_pay(I_pay)  - cons_floor_at(jage, jmar,je,jh)   
  
  !find lower bound 
          jsav_upper=0
          do while ( MAX_SAV .ge. sav_sim(jsav_upper+1) )
           jsav_upper=jsav_upper +1
           if (jsav_upper.EQ.n_sav_sim) exit
          end do 
          
               if (jsav_upper.lt.1) then
                 PRINT *,'Error'
               end if 
                    

               
  !VALUE ASSOCIATED WITH CONSUMING on the consumption floor and saving max amount: V_cf_trpay
         call func_uti(jage,jtrpay,je,jh,jo,jtype4,jtype2,cons_floor(jage, jmar,je,jh),UTI) 
         V_death =  Bequest_val(MAX_SAV+jassets)  + cost_death
         !Calculate expected value associated with MAX_SAV assets tomorrow
         EXPECTED_B_PRIME=0.0d0
         
         if (jage.eq.(n_ret-1)) then 

                    RET_PRIME_INT(:,:,:,:,:)=0.0d0
                    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1, n_MEcat; do jmar_p=1, n_mar
                    !interpolate for assets
                        call InterpolateScalar(as(1:n_as,jtype1,jage+1), RET_PRIME(jhc_ret,jMEcat_p,1:n_as,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p), (MAX_SAV+jassets), RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) )
                     EXPECTED_B_PRIME = EXPECTED_B_PRIME +Tp_h(jh_p, jh,je,jht,jage, jtype2,I_treat,jmc)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p)* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p)* p_surv(jh_p,je,jmar_p, jage+1) + V_death * (1.0-p_surv(jh_p, je,jmar_p, jage+1)) )  
                    end do; end do; end do  ; end do; end do
  
         else 
             EXP_L_prime_INT(:,:, :,:,:)=0.0d0
                    do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar;  do jzhc_p=1,n_zhc
                    call TYPE3_DET(jh_p,jr_p,jtype3_p) 
      !interpolate with respect to human capital and assets    
                     if (jo_act.eq.1) then
                    call interpol2d(as(1:n_as,jtype1,jage+1), xfn(   1: n_hc_grid , jage+1), EXP_L_prime(1:n_as,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,1), (MAX_SAV+jassets), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p))                           
                     else 
                    call interpol2d(as(1:n_as,jtype1,jage+1), xfn(   1: n_hc_grid , jage+1), EXP_L_prime(1:n_as,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,2), (MAX_SAV+jassets), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p))                           

                                             end if 
                    
                    if ((jtrpay.eq.1).OR.(jtrpay.eq.2)) then ! treat 
                    EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(1,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))
     else 
                             EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(jmc+1,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))

     end if 
     
     
     end do ; end do ; end do   ; end do; end do  
  
         end if 
                 
         V_cf_trpay = UTI + beta_hat(je)* EXPECTED_B_PRIME  

         
 !VALUE ASSOCIATED WITH CONSUMING ALL AND 0 ASSETS: V_c_all_trpay
         
           call func_uti(jage,jtrpay,je,jh,jo,jtype4,jtype2,(base_1_pay(I_pay)/(1.0d0+taxc)),UTI) 
         V_death =  Bequest_val(0.0d0)  + cost_death
         !Calculate expected value associated with MAX_SAV assets tomorrow
         EXPECTED_B_PRIME=0.0d0
         
         if (jage.eq.(n_ret-1)) then 

                    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1, n_MEcat; do jmar_p=1, n_mar
                     EXPECTED_B_PRIME = EXPECTED_B_PRIME + Tp_h(jh_p, jh,je,jht,jage, jtype2,I_treat,jmc)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p)* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (RET_PRIME(jhc_ret,jMEcat_p,1,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p) * p_surv(jh_p,je,jmar_p, jage+1) +  V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)) )                                      
                    end do; end do; end do  ; end do; end do
  
         else 
             EXP_L_prime_INT(:,:, :,:,:)=0.0d0
                    do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar;  do jzhc_p=1,n_zhc
                    call TYPE3_DET(jh_p,jr_p,jtype3_p) 
      !interpolate with respect to human capital. Assets will be 0 next period.  
                    if (jo_act.eq.1) then
     call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,1), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )
                    else
     call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,2), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )
                    end if 
                    
                        
          if ((jtrpay.eq.1).OR.(jtrpay.eq.2)) then ! treat 
     EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(1,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))
          else 
                   EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(jmc+1,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))

          end if 
          
          
          end do ; end do ; end do   ; end do; end do  
  
          end if 
                 
         V_c_all_trpay = UTI + beta_hat(je)* EXPECTED_B_PRIME 
         
          
!Iterate over possible savings and find optimal   
         maxB =V_cf_trpay
         jas_p_act=MAX_SAV+jassets
         jcons_act=cons_floor(jage, jmar,je,jh)
do jsav= jsav_upper , 1, -1
! get ASP, CON, and UTI associated with this jsav    
        ASP=  jassets+  sav_sim(jsav)
     if (ASP.le.0.d0)  exit   ! I already did the case where you consume all and save 0, so that's why we have "le," not "lt."
        V_death= cost_death + Bequest_val(ASP)  
        CON= ( income_bt + inc_spouse(jage,je,jh,jtype4,jfix) - base_2_pay(I_pay) -sav_sim(jsav) )/(1.0d0+taxc)
        call  func_uti(jage,jtrpay,je,jh,jo_act,jtype4,jtype2,CON,UTI)           
 
!Solve continuation value given this jsav and jtrpay
          EXPECTED_B_PRIME=0.0d0
         if (jage.eq.(n_ret-1)) then !need to interpolate only for assets
                    RET_PRIME_INT(:,:,:,:,:)=0.0d0


                    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1, n_MEcat; do jmar_p=1, n_mar
                        call InterpolateScalar(as(1:n_as,jtype1,jage+1), RET_PRIME(jhc_ret,jMEcat_p,1:n_as,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p), ASP, RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) )
                     EXPECTED_B_PRIME = EXPECTED_B_PRIME +Tp_h(jh_p, jh,je,jht,jage, jtype2,I_treat,jmc)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p)* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p)* p_surv(jh_p,je,jmar_p, jage+1) +  V_death * (1.0-p_surv(jh_p, je,jmar_p, jage+1)) )                                      
                    end do; end do; end do  ; end do; end do
  
         else 
             EXP_L_prime_INT(:,:, :,:,:)=0.0d0
                    do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar;  do jzhc_p=1,n_zhc
                    call TYPE3_DET(jh_p,jr_p,jtype3_p) 
                    if (jo_act.eq.1) then
                    call interpol2d(as(1:n_as,jtype1,jage+1), xfn(   1: n_hc_grid , jage+1), EXP_L_prime(1:n_as,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,1), ASP, HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p))                           
                    else 
                    call interpol2d(as(1:n_as,jtype1,jage+1), xfn(   1: n_hc_grid , jage+1), EXP_L_prime(1:n_as,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,2), ASP, HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p))                           
                    end if 
                    
                    
                    
                    if ((jtrpay.eq.1).OR.(jtrpay.eq.2)) then ! treat 
                    EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(1,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act)  *(EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p)* p_surv(jh_p,je,jmar_p, jage+1) + V_death * (1.0-p_surv(jh_p,je,jmar_p, jage+1)) ) 
                    else 
                    EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(jmc+1,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act)  *(EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p)* p_surv(jh_p,je,jmar_p, jage+1) + V_death * (1.0-p_surv(jh_p,je,jmar_p, jage+1)) ) 

                    end if 
                    
                    
                    end do ; end do ; end do   ; end do; end do  
  
          end if 
                                 
B_act=UTI+ beta_hat(je)*EXPECTED_B_PRIME 


     if (B_act.lt.maxB) exit
              maxB=B_act  
              jas_p_act=ASP
              jcons_act=CON
    
       
 end do !jsav
! Found jcons_act , jas_p_act and B_act assocaited with optimal savings. 
   
! now also need to compare with:    V_c_all_trpay. And if this is max, then set jas and jcons accordingly.
  if  (V_c_all_trpay.gt.maxB) then
      maxB =V_c_all_trpay
      jas_p_act=0.0d0
      jcons_act= base_1_pay(I_pay)/(1.0d0+taxc)
  end if 
  

end if ! for being on the consumption floor

! save CONS, ASP, TR and B for the 3 jtrpay options. 

  CON_sim(jtrpay)=jcons_act
  ASP_sim(jtrpay)=jas_p_act
  TR_sim(jtrpay)=jtr_act
  B_sim(jtrpay)=maxB
  
    end do ! jtrpay
    
    
    ! determine best option and reset I_treat and I_pay
if ((B_sim(1).eq.0.0d0).or.(B_sim(3).eq.0.0d0)) then
PRINT *,'Error'
pause
end if 
if (B_sim(1).ge.B_sim(3)) then
    jtrpay = 1 
    else
    jtrpay = 3 
end if 

jtr_act= TR_sim(jtrpay)
jas_p_act=ASP_sim(jtrpay)
jcons_act=CON_sim(jtrpay)


 else if ((joption_pay.eq.1).AND.(juntreated.EQ.1)) then ! you cannot treat

    
    jtrpay=3 ! you don't treat
    
! I_treat and I_pay    
    
!jtrpay=1 	pay and get treated
!jtrpay=2	don't pay, get treated
!jtrpay=3	don't pay, not treated
    

        I_treat=1 ! not treat
        I_pay=2 ! not pay


!initiate
  CON=0.0d0
  ASP=0.d0
  UTI=0.d0
  TR=0.0d0
  MAX_SAV=0.0d0 
        jcons_act=0.0d0
        !jsav_act=0.0d0 
        jas_p_act=0.0d0
        jtr_act=0.0d0
        B_act=0.0d0
        
!If you are on the consumption floor        
if ( base_1_pay(I_pay) .le. (cons_floor_at(jage, jmar,je,jh))) then ! you get transfer. you only get to this point if you are not on transfer if you don't pay. so you'll still have to calculate the value if you don't pay
 jcons_act=cons_floor(jage, jmar,je,jh)  ! cons floor is at the HH level; CON is after tax HH consumption - it's what goes into utility fn
 jas_p_act=0.d0 
 jtr_act=  cons_floor_at(jage, jmar,je,jh) - base_1_pay(I_pay) 
 call  func_uti(jage,jtrpay,je,jh,jo_act,jtype4,jtype2,cons_floor(jage, jmar,je,jh),UTI)
 V_death =  Bequest_val(0.0d0)  + cost_death 
 
!continuation value
 EXPECTED_B_PRIME=0.0d0
         if (jage.eq.(n_ret-1)) then 


                    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1, n_MEcat; do jmar_p=1, n_mar
                     EXPECTED_B_PRIME = EXPECTED_B_PRIME + Tp_h(jh_p, jh,je,jht,jage, jtype2,I_treat,jmc)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p)* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (RET_PRIME(jhc_ret,jMEcat_p,1,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p) * p_surv(jh_p,je,jmar_p, jage+1) +  V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)) )                                      
                    end do; end do; end do  ; end do; end do
  
         else 
             EXP_L_prime_INT(:,:, :,:,:)=0.0d0
                    do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar;  do jzhc_p=1,n_zhc
                    call TYPE3_DET(jh_p,jr_p,jtype3_p) 
      !interpolate with respect to human capital. Assets will be 0 next period.   
                    if (jo_act.eq.1) then
     call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,1), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )
                    else 
     call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,2), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )

                    end if 
                    
     if ((jtrpay.eq.1).OR.(jtrpay.eq.2)) then ! treat 
     EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(1,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))
     else 
              EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(jmc+1,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))
     end if 
     
     
     end do ; end do ; end do   ; end do; end do  
  
          end if 
                                 
maxB=UTI+ beta_hat(je)*EXPECTED_B_PRIME 
 
else
  jtr_act=0.d0
! find maximum savings on the grid (jsav_upper)
  MAX_SAV= income_bt + inc_spouse(jage,je,jh,jtype4,jfix) -base_2_pay(I_pay)  - cons_floor_at(jage, jmar,je,jh)   
  
  !find lower bound 
          jsav_upper=0
          do while ( MAX_SAV .ge. sav_sim(jsav_upper+1) )
           jsav_upper=jsav_upper +1
           if (jsav_upper.EQ.n_sav_sim) exit
          end do 
          
               if (jsav_upper.lt.1) then
                 PRINT *,'Error'
               end if 
                    

               
  !VALUE ASSOCIATED WITH CONSUMING on the consumption floor and saving max amount: V_cf_trpay
         call func_uti(jage,jtrpay,je,jh,jo,jtype4,jtype2,cons_floor(jage, jmar,je,jh),UTI) 
         V_death =  Bequest_val(MAX_SAV+jassets)  + cost_death
         !Calculate expected value associated with MAX_SAV assets tomorrow
         EXPECTED_B_PRIME=0.0d0
         
         if (jage.eq.(n_ret-1)) then 

                    RET_PRIME_INT(:,:,:,:,:)=0.0d0
                    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1, n_MEcat; do jmar_p=1, n_mar
                    !interpolate for assets
                        call InterpolateScalar(as(1:n_as,jtype1,jage+1), RET_PRIME(jhc_ret,jMEcat_p,1:n_as,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p), (MAX_SAV+jassets), RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) )
                     EXPECTED_B_PRIME = EXPECTED_B_PRIME +Tp_h(jh_p, jh,je,jht,jage, jtype2,I_treat,jmc)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p)* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p)* p_surv(jh_p,je,jmar_p, jage+1) + V_death * (1.0-p_surv(jh_p, je,jmar_p, jage+1)) )  
                    end do; end do; end do  ; end do; end do
  
         else 
             EXP_L_prime_INT(:,:, :,:,:)=0.0d0
                    do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar;  do jzhc_p=1,n_zhc
                    call TYPE3_DET(jh_p,jr_p,jtype3_p) 
      !interpolate with respect to human capital and assets 
                    if (jo_act.eq.1) then
                    call interpol2d(as(1:n_as,jtype1,jage+1), xfn(   1: n_hc_grid , jage+1), EXP_L_prime(1:n_as,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,1), (MAX_SAV+jassets), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p))                           
                    else 
                    call interpol2d(as(1:n_as,jtype1,jage+1), xfn(   1: n_hc_grid , jage+1), EXP_L_prime(1:n_as,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,2), (MAX_SAV+jassets), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p))                           
                    end if 
                    
                    
                    if ((jtrpay.eq.1).OR.(jtrpay.eq.2)) then ! treat 
                    EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(1,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))
     else 
                             EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(jmc+1,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))

     end if 
     
     
     end do ; end do ; end do   ; end do; end do  
  
         end if 
                 
         V_cf_trpay = UTI + beta_hat(je)* EXPECTED_B_PRIME  

         
 !VALUE ASSOCIATED WITH CONSUMING ALL AND 0 ASSETS: V_c_all_trpay
         
           call func_uti(jage,jtrpay,je,jh,jo,jtype4,jtype2,(base_1_pay(I_pay)/(1.0d0+taxc)),UTI) 
         V_death =  Bequest_val(0.0d0)  + cost_death
         !Calculate expected value associated with MAX_SAV assets tomorrow
         EXPECTED_B_PRIME=0.0d0
         
         if (jage.eq.(n_ret-1)) then 

                    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1, n_MEcat; do jmar_p=1, n_mar
                     EXPECTED_B_PRIME = EXPECTED_B_PRIME + Tp_h(jh_p, jh,je,jht,jage, jtype2,I_treat,jmc)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p)* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (RET_PRIME(jhc_ret,jMEcat_p,1,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p) * p_surv(jh_p,je,jmar_p, jage+1) +  V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)) )                                      
                    end do; end do; end do  ; end do; end do
  
         else 
             EXP_L_prime_INT(:,:, :,:,:)=0.0d0
                    do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar;  do jzhc_p=1,n_zhc
                    call TYPE3_DET(jh_p,jr_p,jtype3_p) 
      !interpolate with respect to human capital. Assets will be 0 next period.     
                    if (jo_act.eq.1) then
     call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,1), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )
                    else 
                             call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,2), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )

                    end if 
                    
          if ((jtrpay.eq.1).OR.(jtrpay.eq.2)) then ! treat 
     EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(1,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))
          else 
                   EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(jmc+1,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))

          end if 
          
          
          end do ; end do ; end do   ; end do; end do  
  
          end if 
                 
         V_c_all_trpay = UTI + beta_hat(je)* EXPECTED_B_PRIME 
         
          
!Iterate over possible savings and find optimal   
         maxB =V_cf_trpay
         jas_p_act=MAX_SAV+jassets
         jcons_act=cons_floor(jage, jmar,je,jh)
do jsav= jsav_upper , 1, -1
! get ASP, CON, and UTI associated with this jsav    
        ASP=  jassets+  sav_sim(jsav)
     if (ASP.le.0.d0)  exit   ! I already did the case where you consume all and save 0, so that's why we have "le," not "lt."
        V_death= cost_death + Bequest_val(ASP)  
        CON= ( income_bt + inc_spouse(jage,je,jh,jtype4,jfix) - base_2_pay(I_pay) -sav_sim(jsav) )/(1.0d0+taxc)
        call  func_uti(jage,jtrpay,je,jh,jo_act,jtype4,jtype2,CON,UTI)           
 
!Solve continuation value given this jsav and jtrpay
          EXPECTED_B_PRIME=0.0d0
         if (jage.eq.(n_ret-1)) then !need to interpolate only for assets; for human capital, we discretize it into bigger categories, so no need to interpolate.
                    RET_PRIME_INT(:,:,:,:,:)=0.0d0


                    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1, n_MEcat; do jmar_p=1, n_mar
                        call InterpolateScalar(as(1:n_as,jtype1,jage+1), RET_PRIME(jhc_ret,jMEcat_p,1:n_as,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p), ASP, RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) )
                     EXPECTED_B_PRIME = EXPECTED_B_PRIME +Tp_h(jh_p, jh,je,jht,jage, jtype2,I_treat,jmc)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p)* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p)* p_surv(jh_p,je,jmar_p, jage+1) +  V_death * (1.0-p_surv(jh_p, je,jmar_p, jage+1)) )                                      
                    end do; end do; end do  ; end do; end do
  
         else 
             EXP_L_prime_INT(:,:, :,:,:)=0.0d0
                    do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar;  do jzhc_p=1,n_zhc
                    call TYPE3_DET(jh_p,jr_p,jtype3_p) 
                    if (jo_act.eq.1) then
                    call interpol2d(as(1:n_as,jtype1,jage+1), xfn(   1: n_hc_grid , jage+1), EXP_L_prime(1:n_as,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,1), ASP, HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p))                           
                    else 
                    call interpol2d(as(1:n_as,jtype1,jage+1), xfn(   1: n_hc_grid , jage+1), EXP_L_prime(1:n_as,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,2), ASP, HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p))                           

                    end if 
                    
                                            
                    if ((jtrpay.eq.1).OR.(jtrpay.eq.2)) then ! treat 
                    EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(1,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act)  *(EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p)* p_surv(jh_p,je,jmar_p, jage+1) + V_death * (1.0-p_surv(jh_p,je,jmar_p, jage+1)) ) 
                    else 
                    EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(jmc+1,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act)  *(EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p)* p_surv(jh_p,je,jmar_p, jage+1) + V_death * (1.0-p_surv(jh_p,je,jmar_p, jage+1)) ) 

                    end if 
                    
                    
                    end do ; end do ; end do   ; end do; end do  
  
          end if 
                                 
B_act=UTI+ beta_hat(je)*EXPECTED_B_PRIME 


     if (B_act.lt.maxB) exit
              maxB=B_act  
              jas_p_act=ASP
              jcons_act=CON
    
       
 end do !jsav
! Found jcons_act , jas_p_act and B_act assocaited with optimal savings. 
   
! now also need to compare with:    V_c_all_trpay. And if this is max, then set jas and jcons accordingly.
  if  (V_c_all_trpay.gt.maxB) then
      maxB =V_c_all_trpay
      jas_p_act=0.0d0
      jcons_act= base_1_pay(I_pay)/(1.0d0+taxc)
  end if 
  

end if ! for being on the consumption floor

! save CONS, ASP, TR and B for the 3 jtrpay options. 

  CON_sim(jtrpay)=jcons_act
  ASP_sim(jtrpay)=jas_p_act
  TR_sim(jtrpay)=jtr_act
  B_sim(jtrpay)=maxB
  
  
  
jtr_act= TR_sim(jtrpay)
jas_p_act=ASP_sim(jtrpay)
jcons_act=CON_sim(jtrpay) 
    
    

    
    
end if ! for having the options in trpay


    if (jtrpay.eq.1) then
        I_treat=2 ! treat
        I_pay=1 ! pay
    else if (jtrpay.eq.2) then
        I_treat=2 ! treat
        I_pay=2 ! not pay
    else if (jtrpay.eq.3) then
        I_treat=1 ! not treat
        I_pay=2 ! not pay
    end if 

    

    
end if ! for being on consumption floor even if you don't pay
end if ! for having a health shock

!If you don't get treated, then there is a probability the shock and R don't get reported

!initialize to not having shock
jdp_new=jdp
jdu_new=jdu
js_new=js
jr_rep=jr

!dp (2,4,6,8)
!if ((jtype2.eq.2).OR.(jtype2.eq.4).OR.(jtype2.eq.6).OR.(jtype2.eq.8) ) then
if (jdp.eq.2) then
    rb=ran0(idum) ! draw a random # between 0-1
    if (rb.gt.neu_shocks(jh,1,I_treat))   then
      jdp_new=1 
    end if 
end if 


!du (5,6,7,8)
!if ((jtype2.eq.5).OR.(jtype2.eq.6).OR.(jtype2.eq.7).OR.(jtype2.eq.8) ) then
 if (jdu.eq.2) then
    rb=ran0(idum) ! draw a random # between 0-1
    if (rb.gt.neu_shocks(jh,2,I_treat))   then
      jdu_new=1 
    end if 
end if

!s (3,4,7,8)
 if (js.eq.2) then
    rb=ran0(idum) ! draw a random # between 0-1
    if (rb.gt.neu_shocks(jh,3,I_treat))   then
      js_new=1 
    end if 
 end if
 
 ! R
 if (jr.eq.2) then
         rb=ran0(idum) ! draw a random # between 0-1
    if (rb.gt.neu_R(jo_act))   then
      jr_rep=1 
    end if 
 end if 
 
 ! but if you get a gov transfer (have Medicaid), then all shocks and R reported correctly; this should already be OK for shocks. making sure ok for R too.
 if (jtr_act.gt.0.0d0) then 
 jdp_new=jdp
 jdu_new=jdu
 js_new=js
 jr_rep=jr
 end if 
 

call TYPE2_DET(jdu_new,jdp_new,js_new,jtype2_new)

!Write to file
 write(12,1008) jsample, jage, je, jfix, jht, jmar, jh,  jr, jtype2, jtype4, jo ,   jo_act,  jzt, jx_int, jx_int_p, oop(jtype2, jh, jage,jMEcat,jo_act)*5200.0d0 , Charges(jtype2, jh, jage,jMEcat)*5200.0d0, jcons_act*5200.0d0, jassets*5200.0d0 , jas_p_act*5200.0d0,  jtr_act*5200.0d0  , wage_penalty(jo,je)*Z_INT_OFFER   , 5200.0d0*actual_earn, max(0.0d0, (hrs_offer(jo_act) -phi(je,jh,jtype2) ))*100.0d0, exp( log(wage_penalty(jo_act,je)*Z_INT_OFFER)+random_normal()*stdze(je) ) ,  inc_spouse(jage,je,jh,jtype4,jfix)*5200.0d0, jtrpay  , oop_spouse(jage,je,jtype4,jo_act)*5200.0d0 ,base_1_pay(1)*5200, base_1_pay(2)*5200, base_2_pay(1)*5200, base_2_pay(2)*5200, 0, jtype2_new, jr_rep, joption_pay, juntreated, B_sim(jtrpay)
 
 
 
!draw marital status for next period 
            rf=ran0(idum) ! draw a random # between 0-1
            
                        if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.0.AND.jage.LE.40) then
            write(14,1005) rf 
            else if (jage.LE.40) then
 
        rf=ran_hypo(12+11*(jage-1))
                                end if
                                
            if (rf.lt.(Tp_mar(1,jmar,je,jage,jh,jinc,jo_act))) then
            jmar_p=1
            else 
            jmar_p=2
            end if
          jmar=jmar_p    
!Health Risk R
          rb=ran0(idum) ! draw a random # between 0-1
                                  if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.0.AND.jage.LE.40) then
            write(14,1005) rb 
            else if (jage.LE.40) then 
        rb=ran_hypo(13+11*(jage-1))
                                end if
          if (rb.lt.p_risk(jr, 1,jage,jh)) then
            jr_p=1
            else 
            jr_p=2
            end if
          jr=jr_p  
          
          
                       if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.6) then 
                       jr=1
                       end if
                       
              if (rb.lt.p_risk(jr_alt, 1,jage,jh_alt)) then
            jr_p_alt=1
            else 
            jr_p_alt=2
            end if
          jr_alt=jr_p_alt                    
                       
                       
!Health
          rb=ran0(idum) ! draw a random # between 0-1
                                  if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.0.AND.jage.LE.40) then
            write(14,1005) rb 
            else if (jage.LE.40) then 
        rb=ran_hypo(14+11*(jage-1))
                                end if
                                
if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.9.AND.jage.EQ.6 ) then ! age 30, do not allow H to go down. just keep H the same
	   jh=jh
else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.10.AND.jage.EQ.16 ) then ! age 40, do not allow for d shocks; just keep H the same
  	   jh=jh
else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.11.AND.jage.EQ.26 ) then ! age 50, do not allow for d shocks; just keep H the same
 	   jh=jh
else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.12.AND.jage.EQ.36 ) then ! age 60, do not allow for d shocks; just keep H the same
 	   jh=jh
    
else if (jExperimentNumber.EQ.0.AND.(jBenchmark_exp.EQ.13.OR.jBenchmark_exp.EQ.17.OR.jBenchmark_exp.EQ.29).AND.jage.EQ.6.AND.(jh.EQ.2.or.jh.EQ.3) ) then ! age 30, ALL GET du shocks and H goes down by 1
       if (jh.eq.3) then
       jh=2
       else if (jh.eq.2) then
       jh=1 
       end if       
else if (jExperimentNumber.EQ.0.AND.(jBenchmark_exp.EQ.14.OR.jBenchmark_exp.EQ.18.OR.jBenchmark_exp.EQ.27.OR.jBenchmark_exp.EQ.30).AND.jage.EQ.16 ) then ! age 40, ALL GET du shocks and H goes down by 1
        if (jh.eq.3) then
       jh=2
       else if (jh.eq.2) then
       jh=1 
       end if  
else if (jExperimentNumber.EQ.0.AND.(jBenchmark_exp.EQ.15.OR.jBenchmark_exp.EQ.19.OR.jBenchmark_exp.EQ.31).AND.jage.EQ.26 ) then ! age 50, ALL GET du shocks and H goes down by 1
        if (jh.eq.3) then
       jh=2
       else if (jh.eq.2) then
       jh=1 
       end if  
else if (jExperimentNumber.EQ.0.AND.(jBenchmark_exp.EQ.16.OR.jBenchmark_exp.EQ.20.OR.jBenchmark_exp.EQ.32).AND.jage.EQ.36 ) then ! age 60, ALL GET du shocks and H goes down by 1
           if (jh.eq.3) then
       jh=2
       else if (jh.eq.2) then
       jh=1 
       end if  
      
else                                 
                                
                                
            bb=1
            prob=Tp_h(bb, jh,je,jht,jage, jtype2,I_treat,jmc) 
         92 if (rb.lt.prob) then
            jh_p=bb
            else 
            bb=bb+1
                if (bb.le.n_h) then
                prob=prob+Tp_h(bb, jh,je,jht,jage, jtype2,I_treat,jmc)
                go to 92
                else 
              jh_p= n_h  
                end if
            end if  
            jh=jh_p           
end if 
           !Health if no health shocks - set jtype2 to 1. 
            bb=1
            prob=Tp_h(bb, jh_alt,je,jht,jage, 1,I_treat,jmc) 
        292 if (rb.lt.prob) then
            jh_alt=bb
            else 
            bb=bb+1
                if (bb.le.n_h) then
                prob=prob+Tp_h(bb, jh_alt,je,jht,jage, 1,I_treat,jmc)
                go to 292
                else 
              jh_alt= n_h  
                end if
            end if 
            
            

!Next period states needed only for those of working age, but not 64: 
if (jage.NE.(n_ret-1)) then    
             
        jx_int=max(0.0d0, min(jx_int_p , dble(n_hc)) )
        

       
!Draw employment offer jo

      ro=ran0(idum) ! draw a random # between 0-1
                        if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.0.AND.jage.LE.39) then
            write(14,1005) ro 
            else if (jage.LE.39) then
        ro=ran_hypo(15+11*(jage-1))
                                      end if
                                      
        bb=1
        prob=Pr_o_i(bb,jo_act, je,jage)
 
     34 if (ro.lt.prob) then
        jo=bb
        else 
        bb=bb+1
        if (bb.le.n_o) then
        prob=prob+Pr_o_i(bb,jo_act, je,jage)

            
        go to 34
        else 
            jo=n_o
        end if
        end if   
        

        
        
!Draw zt next period
                               rb=ran0(idum) ! draw a random # between 0-1
                                                       if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.0.AND.jage.LE.39) then
            write(14,1005) rb 
            else if (jage.LE.39) then 

        rb=ran_hypo(16+11*(jage-1))
                                      end if
                                      
        if (jo_act.eq.1) then
         
		            bb=1
		            prob=p_zt(je,bb,1)
89		            if (rb.lt.prob) then
		            jzt=bb
		            else 
		            bb=bb+1
		            if (bb.le.n_zt) then
		            prob=prob+p_zt(je,bb,1)
		            go to 89
                    else 
                        jzt=n_zt
		            end if
                    end if  
                    
 
        else 

		            bb=1
		            prob=p_zt(je,bb,2)
99		            if (rb.lt.prob) then
		            jzt=bb
		            else 
		            bb=bb+1
		            if (bb.le.n_zt) then
		            prob=prob+p_zt(je,bb,2)
		            go to 99
                    else 
                        jzt=n_zt
		            end if
                    end if  
        end if 
        
!Spouse's work status next period (jos)

                        ! draw spouse work status  1=not work, 2=work ; if not married, doesn't matter, jos will be drawn but not used.
                        rb=ran0(idum) 
                                                                               if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.0.AND.jage.LE.39) then
            write(14,1005) rb 
            else if (jage.LE.39) then

        rb=ran_hypo(17+11*(jage-1))
                                      end if
                        if (rb.lt.(prob_os(jage+1,je,jh,1,jfix))) then ! this is the jh next period   
                        jos=1
                        else 
                        jos=2
                        end if 
               
end if ! for age not being 64
end if ! for jage<n_ret       


		
!****************************************************************************************************************************************


if (jage.ge.n_ret) then !individual is retired	
           
!Draw health shocks jtype2
                  rf=ran0(idum) ! draw a random # between 0-1
		            bb=1
		            prob=shocks_risk(je,bb, jh, jr, jage) 
		         37 if (rf.lt.prob) then
		            jtype2=bb
		            else 
		            bb=bb+1
		            if (bb.le.n_type2) then
		            prob=prob+shocks_risk(je,bb, jh, jr, jage)
		            go to 37
                    else 
                        jtype2=n_type2
		            end if
                    end if
                    
                    
                     if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.1 ) then
 if (jtype2.eq.3.OR.jtype2.eq.4.OR.jtype2.eq.7.OR.jtype2.eq.8) then 
     jtype2=jtype2-2
end if
else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.2 ) then
 if (jtype2.eq.5.OR.jtype2.eq.6.OR.jtype2.eq.7.OR.jtype2.eq.8) then 
     jtype2=jtype2-4
end if	            
else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.3 ) then
 if (jtype2.eq.2.OR.jtype2.eq.4.OR.jtype2.eq.6.OR.jtype2.eq.8) then
     jtype2=jtype2-1
end if		         
else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.4 ) then
 if (jtype2.eq.3.OR.jtype2.eq.4)   then 
     jtype2=jtype2-2
 else if (jtype2.eq.5.OR.jtype2.eq.6) then 
      jtype2=jtype2-4
else if (jtype2.eq.7.OR.jtype2.eq.8) then 
       jtype2=jtype2-6
end if 
else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.5 ) then
 jtype2=1
end if	


!Draw survival		            
                     rx=ran0(idum)
		             if (rx.gt.p_surv(jh,je,jmar , jage)) then !individual dies	
                     go to 33 !simulate new individual 
                     end if	             
!Draw the medical expenditure catastrophic shock
                  rf=ran0(idum) ! draw a random # between 0-1
		            bb=1
		            prob= Prob_MEcat(bb)
		         18 if (rf.lt.prob) then
		            jMEcat=bb
		            else 
		            bb=bb+1
		            if (bb.le.n_MEcat) then
		            prob=prob+Prob_MEcat(bb)
		            go to 18
                    else 
                        jMEcat=n_MEcat
		            end if
                    end if            
                     
		             


!income
income = max(0.0d0, (interest-1.0d0)*jassets + pen(jhc_ret,  je,jfix) + inc_spouse(jage,je,1,jmar,jfix) -  prem_mcare(jmar)  - max(0.d0,  oop(jtype2, jh, jage,jMEcat,1) +  oop_spouse(jage,je,jmar,1) - 0.075d0*((interest-1.0d0)*jassets + pen(jhc_ret,  je,jfix) + inc_spouse(jage,je,1,jmar,jfix) ) )) !taxable income. the last term is OOP in excess of a threshold of your income, which is deductable.
base_2 = prem_mcare(jmar) + oop(jtype2, jh, jage,jMEcat,1) + oop_spouse(jage,je,jmar,1)  + max(0.0d0,  (( (income*5200.0d0) - tax_lambda(jmar) *( (income*5200.0d0)**(1.0d0-tax_tau(jmar))) )/5200 )) ! total resources subtracted: premiums, OOP, and tax, at the HH level
                
        
!Initialize
 RET_act=0.0d0
 jtr_act=0.0d0
 jcons_act=0.0d0
 jas_p_act=0.0d0
 maxB=0.0d0
  
if (((interest*jassets+ pen(jhc_ret,  je,jfix)  + inc_spouse(jage,je,1,jmar,jfix) - base_2 )) .le. (cons_floor_at(jage, jmar,je,jh)) ) then    ! get transfer                                                                                                                                  
    jtr_act= cons_floor_at(jage, jmar,je,jh) - (interest*jassets+ pen(jhc_ret, je,jfix) + inc_spouse(jage,je,1,jmar,jfix) -base_2) 
    jcons_act=cons_floor(jage, jmar,je,jh)
    jas_p_act=0.0d0

else
   jtr_act=0.d0
!Find max savings point
   MAX_SAV=  interest*jassets  + pen(jhc_ret, je,jfix) + inc_spouse(jage,je,1,jmar,jfix) -base_2  -cons_floor_at(jage, jmar,je,jh) -jassets
   
!Find lower bound
            jsav_upper_ret=0
          do while ( MAX_SAV .ge. sav_sim(jsav_upper_ret+1) )
           jsav_upper_ret=jsav_upper_ret +1
           if (jsav_upper_ret.EQ.n_sav_sim) exit
          end do 
          
               if (jsav_upper_ret.lt.1) then
                 PRINT *,'Error'
               end if 
               
          
!Find value of consuming all and saving nothing
   CON=((interest*jassets+ pen(jhc_ret,  je,jfix)  + inc_spouse(jage,je,1,jmar,jfix) - base_2 )) /(1.0d0+taxc)      
        call func_uti_ret(je,jh,jmar,jtype2,CON,UTI)
                        V_death =  Bequest_val(0.0d0)  + cost_death
         !Calculate expected value associated with MAX_SAV assets tomorrow
         
                if (jage.eq.n_age) then
                 V_c_all_trpay=UTI+ beta_hat(je)* V_death  
                else 
                   EXPECTED_RET_PRIME=0.0d0
                    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1, n_MEcat; do jmar_p=1, n_mar
                     EXPECTED_B_PRIME = EXPECTED_B_PRIME + Tp_h(jh_p, jh,je,jht,jage, jtype2,2,1)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p)* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (RET_PRIME(jhc_ret,jMEcat_p,1,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p) * p_surv(jh_p,je,jmar_p, jage+1) +  V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)) )                                      
                    end do; end do; end do  ; end do; end do 
                 
                 V_c_all_trpay = UTI + beta_hat(je)* EXPECTED_B_PRIME 
                end if 
         
!Find value of consuming the consumption floor and saving as much as possible 
      call func_uti_ret(je,jh,jmar,jtype2,cons_floor(jage, jmar,je,jh),UTI)
      ASP= MAX_SAV+jassets
      V_death =  Bequest_val(MAX_SAV+jassets)  + cost_death    
                   
                if (jage.eq.n_age) then
                 V_cf_trpay=UTI+ beta_hat(je)* V_death  
                else 
                  EXPECTED_RET_PRIME=0.0d0   
                do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1,n_MEcat; do jmar_p=1,n_mar
                    call InterpolateScalar(as(1:n_as,jtype1,jage+1), RET_PRIME(jhc_ret,jMEcat_p,1:n_as,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p), ASP, RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) )
                EXPECTED_RET_PRIME= EXPECTED_RET_PRIME+ Tp_h(jh_p, jh,je,jht,jage, jtype2,2,1)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p) *  Tp_mar(jmar_p,jmar,je,jage,1,1,1) * ( p_surv(jh_p,je,jmar_p, jage+1)* RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) + V_death* (1.0 -p_surv(jh_p,je,jmar_p, jage+1) ) )
                end do; end do; end do  ; end do; end do
        V_cf_trpay=UTI+ beta_hat(je)*EXPECTED_RET_PRIME
                end if 
                 
               

!loop for optimal savings
         maxB =V_cf_trpay
         jas_p_act=MAX_SAV+jassets
         jcons_act=cons_floor(jage, jmar,je,jh)
   do jsav= jsav_upper_ret , 1, -1
       
!Find CON, ASP and UTI

    ASP=  jassets+  sav_sim(jsav)
                if (ASP.lt.0.0d0)  exit
                V_death=cost_death + Bequest_val(ASP)
                
        CON= ( ((interest-1.0d0))*jassets  + pen(jhc_ret,  je,jfix)+ inc_spouse(jage,je,1,jmar,jfix) -base_2 -sav_sim(jsav) )/(1.0d0+taxc)

        if (CON.lt.cons_floor(jage, jmar,je,jh)) then 
         PRINT *,'Error'
          pause
        end if      
 
     call func_uti_ret(je,jh,jmar,jtype2,CON,UTI)
     
!Find continuation value
            
                EXPECTED_RET_PRIME=0.0d0 
                
                if (jage.eq.n_age) then
                 RET_act=UTI+ beta_hat(je)* V_death  
                else 
                    
                do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1,n_MEcat; do jmar_p=1,n_mar
                    call InterpolateScalar(as(1:n_as,jtype1,jage+1), RET_PRIME(jhc_ret,jMEcat_p,1:n_as,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p), ASP, RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) )
                EXPECTED_RET_PRIME= EXPECTED_RET_PRIME+ Tp_h(jh_p, jh,je,jht,jage, jtype2,2,1)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p) *  Tp_mar(jmar_p,jmar,je,jage,1,1,1) * ( p_surv(jh_p,je,jmar_p, jage+1)* RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) + V_death* (1.0 -p_surv(jh_p,je,jmar_p, jage+1) ) )
                end do; end do; end do  ; end do; end do
        RET_act=UTI+ beta_hat(je)*EXPECTED_RET_PRIME
                end if 
                
!Find max value
                if (RET_act.lt.maxB) exit
                maxB=RET_act 
                jas_p_act=ASP
                jcons_act=CON         
    		    	    
56 end do !for jsav for retired ppl
   
     if  (V_c_all_trpay.gt.maxB) then
      maxB =V_c_all_trpay
      jas_p_act=0.0d0
      jcons_act= ((interest*jassets+ pen(jhc_ret,  je,jfix)  + inc_spouse(jage,je,1,jmar,jfix) - base_2 )) /(1.0d0+taxc)
     end if 
     
end if ! for being on tr
      


   
!Write decisions to file
write(12,1008)  jsample, jage, je, jfix, jht, jmar, jh,  jr, jtype2, 0, 0 ,  0,  0, 0.0d0 , 0.0d0, oop(jtype2, jh, jage,jMEcat,1)*5200.0d0 , 0.0d0 , jcons_act*5200.0d0, jassets*5200.0d0 , jas_p_act*5200.0d0,  jtr_act*5200.0d0  , 0.0d0   , pen(jhc_ret,  je,jfix)*5200.0d0 ,  0.0d0 ,  0.0d0,  inc_spouse(jage,je,1,jmar,jfix)*5200.0d0, 0  , oop_spouse(jage,je,jmar,1)*5200.0d0, ((interest*jassets+  pen(jhc_ret,  je,jfix)  + inc_spouse(jage,je,1,jmar,jfix) - base_2 ))*5200.0d0 , 0.0, base_2*5200.0d0, 0.0, jhc_ret, 0, 0, 0, 0, maxB

if (jage.eq.n_age) then
go to 33 !simulate new individual      
end if                            
      
! draw marital status for next period
            rf=ran0(idum) ! draw a random # between 0-1
            if (rf.lt.(Tp_mar(1,jmar,je,jage,1,1,1))) then
            jmar_p=1
            else 
            jmar_p=2
            end if
          jmar=jmar_p            
!Risk factors
          rb=ran0(idum) ! draw a random # between 0-1
            bb=1
            prob=p_risk(jr, bb,jage,jh)
         39 if (rb.lt.prob) then
            jr_p=bb
            else 
            bb=bb+1
            if (bb.LE.n_r) then
            prob=prob+p_risk(jr, bb,jage,jh)
            go to 39
            else 
            jr_p=n_r  
            end if
            end if
           jr=jr_p        
!Health
          rb=ran0(idum) ! draw a random # between 0-1
            bb=1
            prob=Tp_h(bb, jh,je,jht,jage, jtype2,2,1) ! YOU ALWAYS GET TREATED IF RETIRED
         32 if (rb.lt.prob) then
            jh_p=bb
            else 
            bb=bb+1
            if (bb.le.n_h) then
            prob=prob+Tp_h(bb, jh,je,jht,jage, jtype2,2,1)
            go to 32
            else 
            jh_p= n_h 
            end if
            end if 
    	    jh=jh_p          
end if !for being retired
           
!****************************************************************************************************************************************            
!Assets next period  
  if (jage.lt.n_age)then  
            jassets=jas_p_act          
  end if                 		            
!****************************************************************************************************************************************
  
end do !age

33 end do !sample


end do; end do; end do ! jht, jfix, je

close (12)
close (13)
close (14)
close (15)
close (16)
end do

End Subroutine